feat: multi-branch decay support (B3) — per-branch transition types, intensity cutoff, raw-value plotting

- Extract per-branch forbiddenness codes from ENSDF (not defaulting to 'A')
- BranchConfig dataclass with endpoint, transition type, intensity validation
- Multi-branch BetaSpectrum with per-branch calculators and universal components
- Raw (un-normalized) plotting so total = arithmetic sum of branches
- --intensity-cutoff CLI option to filter negligible branches
- Branch-aware CSV export with per-branch spectra and components
- 30 comprehensive tests covering validation, calculation, export, edge cases

README.md

@@ -149,5 +149,6 @@ For a complete list of implemented and planned features, see [TODO.md](TODO.md).
 
 - **Physics corrections**: phase space, Fermi function, finite size, screening, exchange, radiative (with delta_cut resummation)
 - **Detector response**: Gaussian, Gaussian+tail, Tikhonov models with tabulated support
+- **Multi-branch decay**: intensity-weighted sum of branch spectra, per-branch calculators with branch-specific endpoint, transition type (from ENSDF), and `--intensity-cutoff` CLI filter
 - **Analysis tools**: χ² fitter, C(W) extraction pipeline, Kurie plot analysis
-- **CLI**: `bs_pnpi` with paceENSDF integration, structured logging, CSV metadata headers
+- **CLI**: `bs_pnpi` with paceENSDF integration, structured logging, CSV metadata headers, branch info in dry-run and CSV output

TODO.md

@@ -101,10 +101,21 @@ ______________________________________________________________________
 
 ### B3. Multiple Branch Decays
 
-- [ ] Add support for multiple branch decays
+- [x] ~~Add support for multiple branch decays~~
-  - Calculate total spectral shape weighted by branch intensities
+  - ~~Calculate total spectral shape weighted by branch intensities~~
-  - Handle decays with significant branching to excited states
+  - ~~Handle decays with significant branching to excited states~~
-  - Sum contributions from all branches with proper weighting
+  - ~~Sum contributions from all branches with proper weighting~~
+  - ~~Implemented: `BranchConfig` dataclass, multi-branch `BetaSpectrum` with per-branch calculators~~
+  - ~~Branch data from paceENSDF: intensity, endpoint (Q-value - excitation energy), transition type~~
+  - ~~Intensity normalization to 100% across all branches~~
+  - ~~CSV export: total spectrum, per-branch spectra, per-branch components (PhaseSpace, Radiative), universal components (Fermi, Screening, etc.)~~
+  - ~~Plot: total spectrum + individual branches with legend (transition type, endpoint energy)~~
+  - ~~Debug plot: vertical layout with per-branch panels, universal components panel~~
+  - ~~Energy grid extends to max branch endpoint, per-branch spectra masked beyond their endpoint~~
+  - ~~Single-branch = trivial case of multi-branch (no separate mode)~~
+  - ~~Per-branch transition types extracted from ENSDF (not defaulting to 'A')~~
+  - ~~`--intensity-cutoff` CLI option to filter out negligible branches~~
+  - ~~Branch spectra plotted as raw (un-normalized) contributions so total = sum of branches~~
 
 ### B4. Beta+ Decay Support
 
@@ -207,9 +218,9 @@ ______________________________________________________________________
 
 ## Current Status
 
-**Version:** 0.3.0\
+**Version:** 0.3.1\
-**Implemented:** Phase space, Fermi function, finite size, screening, exchange, radiative corrections (with delta_cut resummation). Detector response module with analytical models (Gaussian, Gaussian+tail, Tikhonov), convolution API, declarative config integration. χ² curve fitting framework (CurveFitter) with confidence intervals, profile likelihood, and correlation analysis. C(W) shape factor extraction pipeline (CWExtractor) with Kurie plot analysis, parametrized fitting, and g_V/g_A extraction. CLI interface (`bs_pnpi`) with paceENSDF integration, structured logging (-v/-vv/-q), --dry-run, --version, --log-file. CSV metadata headers. JSON input with full detector param support. Comprehensive test suite (195 tests). Notebook quality control with nbmake and auto-save plot hooks.
+**Implemented:** Phase space, Fermi function, finite size, screening, exchange, radiative corrections (with delta_cut resummation). Detector response module with analytical models (Gaussian, Gaussian+tail, Tikhonov), convolution API, declarative config integration. χ² curve fitting framework (CurveFitter) with confidence intervals, profile likelihood, and correlation analysis. C(W) shape factor extraction pipeline (CWExtractor) with Kurie plot analysis, parametrized fitting, and g_V/g_A extraction. Multi-branch decay support with per-branch calculators and intensity-weighted sum. CLI interface (`bs_pnpi`) with paceENSDF integration, structured logging (-v/-vv/-q), --dry-run, --version, --log-file. CSV metadata headers with branch info. JSON input with full detector param support. Comprehensive test suite (226 tests). Notebook quality control with nbmake and auto-save plot hooks.
 
-**Completed:** A2 — detector response function and convolution routines. A4 — fitter routine and C(W) extraction pipeline. B2 — input flexibility (JSON input, ENSDF auto-population). B8 — CLI & output improvements (logging, CSV headers, dry-run, version). B8.2 — CLI refinement: removed --mode, --transition-type, --decay-index flags; clarified units (MeV everywhere); enhanced logging (INFO shows components, DEBUG shows internals); CSV headers use element notation (Tc99 -> Ru99). B8.3 — debug verification framework: `./output/` directory for artifacts, development workflow step 4.5 with `-vv` parameter consistency checks across all components. B8.4 — plot output with ID tracking (commit hash + UTC timestamp) and two display modes (normal spectrum plot with nuclear data header; debug 4-panel view with all corrections).
+**Completed:** A2 — detector response function and convolution routines. A4 — fitter routine and C(W) extraction pipeline. B2 — input flexibility (JSON input, ENSDF auto-population). B3 — multi-branch decay support (BranchConfig, per-branch calculators, intensity weighting, branch-aware CSV export and plotting, per-branch transition types from ENSDF, `--intensity-cutoff` CLI option, raw-value plotting where total = sum of branches). B8 — CLI & output improvements (logging, CSV headers, dry-run, version). B8.2 — CLI refinement: removed --mode, --transition-type, --decay-index flags; clarified units (MeV everywhere); enhanced logging (INFO shows components, DEBUG shows internals); CSV headers use element notation (Tc99 -> Ru99). B8.3 — debug verification framework: `./output/` directory for artifacts, development workflow step 4.5 with `-vv` parameter consistency checks across all components. B8.4 — plot output with ID tracking (commit hash + UTC timestamp) and two display modes (normal spectrum plot with nuclear data header; debug 4-panel view with all corrections).
 
 **Next immediate step:** A3 — data processing pipeline for experimental spectra (background subtraction, energy calibration, dead-time correction, pulse pile-up correction).

beta_spectrum/__init__.py

@@ -11,7 +11,7 @@ from .components.exchange import ExchangeCorrection
 from .components.radiative import RadiativeCorrection
 from .components.detector_response import DetectorResponse
 
-from .spectrum import BetaSpectrum, SpectrumConfig, BetaSpectrumAnalyzer
+from .spectrum import BetaSpectrum, SpectrumConfig, BetaSpectrumAnalyzer, BranchConfig
 from .fitter import CurveFitter, FitConfig, FitResult
 from .cw_extractor import CWExtractor, CWExtractionResult, GVAExtractionResult
 from .nuclear_data import (
@@ -44,6 +44,7 @@ __all__ = [
     "BetaSpectrum",
     "SpectrumConfig",
     "BetaSpectrumAnalyzer",
+    "BranchConfig",
     "CurveFitter",
     "FitConfig",
     "FitResult",

beta_spectrum/cli.py

@@ -156,6 +156,15 @@ def _build_parser() -> argparse.ArgumentParser:
         action="store_true",
         help="Validate input and display resolved configuration without computing spectrum.",
     )
+    log_group.add_argument(
+        "--intensity-cutoff",
+        type=float,
+        default=0.0,
+        metavar="FRAC",
+        help="Minimum branch intensity fraction to include in the spectrum "
+        "(default: 0.0 = include all branches). "
+        "Branches below this threshold are excluded from the calculation and plots.",
+    )
 
     return parser
 
@@ -181,10 +190,16 @@ def _run(args: argparse.Namespace) -> None:
     from beta_spectrum.logging_utils import get_logger
     from beta_spectrum.nuclear_data import (
         create_config_from_source,
+        get_decay_info_from_paceENSDF,
         json_to_config,
         load_json_input,
     )
-    from beta_spectrum.spectrum import BetaSpectrum, BetaSpectrumAnalyzer
+    from beta_spectrum.spectrum import (
+        BetaSpectrum,
+        BetaSpectrumAnalyzer,
+        BranchConfig,
+        SpectrumConfig,
+    )
 
     # Determine logging level
     if args.quiet:
@@ -214,16 +229,78 @@ def _run(args: argparse.Namespace) -> None:
     config: SpectrumConfig
     if args.nuclide:
         logger.info("Source: paceENSDF nuclide=%s", args.nuclide)
-        config = create_config_from_source(
+        decay_info = get_decay_info_from_paceENSDF(args.nuclide, "beta_minus")
-            "paceENSDF",
+
-            nuclide=args.nuclide,
+        # Normalize branch intensities to sum to 1.0
-            e_step_MeV=args.e_step,
+        if decay_info.branches:
-            use_detector_response=detector_kwargs.get("use_detector_response", False),
+            total_intensity = sum(b.intensity for b in decay_info.branches)
-            detector_model=detector_kwargs.get("detector_model", "gaussian"),
+            if total_intensity > 0:
-            detector_sigma_a_keV=detector_kwargs.get("detector_sigma_a_keV", 1.0),
+                branches = [
-            detector_tau_keV=detector_kwargs.get("detector_tau_keV", 5.0),
+                    BranchConfig(
-            detector_tail_fraction=detector_kwargs.get("detector_tail_fraction", 0.0),
+                        endpoint_MeV=(
-        )
+                            (decay_info.endpoint_MeV - b.level_energy_keV / 1000.0)
+                            if b.level_energy_keV > 0
+                            else decay_info.endpoint_MeV
+                        ),
+                        transition_type=b.transition_type,
+                        intensity=b.intensity / total_intensity,
+                    )
+                    for b in decay_info.branches
+                    if b.intensity / total_intensity >= args.intensity_cutoff
+                ]
+                config = SpectrumConfig(
+                    Z_parent=decay_info.Z_parent,
+                    Z_daughter=decay_info.Z_daughter,
+                    A_number=decay_info.A_number,
+                    endpoint_MeV=decay_info.endpoint_MeV,
+                    transition_type=decay_info.transition_type,
+                    e_step_MeV=args.e_step,
+                    branches=branches,
+                    intensity_cutoff=args.intensity_cutoff,
+                    use_detector_response=detector_kwargs.get(
+                        "use_detector_response", False
+                    ),
+                    detector_model=detector_kwargs.get("detector_model", "gaussian"),
+                    detector_sigma_a_keV=detector_kwargs.get(
+                        "detector_sigma_a_keV", 1.0
+                    ),
+                    detector_tau_keV=detector_kwargs.get("detector_tau_keV", 5.0),
+                    detector_tail_fraction=detector_kwargs.get(
+                        "detector_tail_fraction", 0.0
+                    ),
+                )
+            else:
+                config = create_config_from_source(
+                    "paceENSDF",
+                    nuclide=args.nuclide,
+                    e_step_MeV=args.e_step,
+                    use_detector_response=detector_kwargs.get(
+                        "use_detector_response", False
+                    ),
+                    detector_model=detector_kwargs.get("detector_model", "gaussian"),
+                    detector_sigma_a_keV=detector_kwargs.get(
+                        "detector_sigma_a_keV", 1.0
+                    ),
+                    detector_tau_keV=detector_kwargs.get("detector_tau_keV", 5.0),
+                    detector_tail_fraction=detector_kwargs.get(
+                        "detector_tail_fraction", 0.0
+                    ),
+                )
+        else:
+            config = create_config_from_source(
+                "paceENSDF",
+                nuclide=args.nuclide,
+                e_step_MeV=args.e_step,
+                use_detector_response=detector_kwargs.get(
+                    "use_detector_response", False
+                ),
+                detector_model=detector_kwargs.get("detector_model", "gaussian"),
+                detector_sigma_a_keV=detector_kwargs.get("detector_sigma_a_keV", 1.0),
+                detector_tau_keV=detector_kwargs.get("detector_tau_keV", 5.0),
+                detector_tail_fraction=detector_kwargs.get(
+                    "detector_tail_fraction", 0.0
+                ),
+            )
         source_type = "paceENSDF"
     elif args.input:
         logger.info("Source: JSON file=%s", args.input)
@@ -241,10 +318,14 @@ def _run(args: argparse.Namespace) -> None:
 
     # Create and calculate spectrum
     start_time = time.time()
+    mode_str = ""
+    if config.branches and len(config.branches) > 1:
+        mode_str = f" (multi-branch: {len(config.branches)} branches)"
     logger.info(
-        "Creating BetaSpectrum for %s, endpoint=%.3f MeV",
+        "Creating BetaSpectrum for %s, endpoint=%.3f MeV%s",
         _decay_notation(config.Z_parent, config.Z_daughter, config.A_number),
         config.endpoint_MeV,
+        mode_str,
     )
 
     spectrum = BetaSpectrum.from_config(config, logger=logger)
@@ -305,6 +386,22 @@ def _print_dry_run_output(config, source_type: str) -> None:
     print(f"Transition: {config.transition_type}")
     print(f"e_step:     {config.e_step_MeV:.4f} MeV")
     print(f"Corrections: {', '.join(enabled) if enabled else 'none'}")
+
+    if config.branches and len(config.branches) > 1:
+        print(f"\nMulti-branch mode: {len(config.branches)} branches")
+        if config.intensity_cutoff > 0:
+            print(
+                f"Intensity cutoff: {config.intensity_cutoff:.4f} ({config.intensity_cutoff*100:.1f}%)"
+            )
+        for i, branch in enumerate(config.branches):
+            print(
+                f"  Branch {i+1}: E₀={branch.endpoint_MeV*1000:.1f} keV, "
+                f"transition={branch.transition_type}, "
+                f"intensity={branch.intensity:.4f}"
+            )
+    else:
+        print("Mode: single-branch")
+
     if config.use_detector_response:
         print(
             f"Detector:   {config.detector_model} (sigma={config.detector_sigma_a_keV} keV, "

beta_spectrum/nuclear_data.py

@@ -94,12 +94,15 @@ class BranchInfo:
         Branch intensity (fraction of decays).
     log_ft : float
         log(ft) value for this branch.
+    transition_type : str
+        Forbiddenness classification for this branch (A, F1, F1U, F2, ...).
     """
 
     level_index: int
     level_energy_keV: float
     intensity: float
     log_ft: Optional[float]
+    transition_type: str = "A"
 
 
 def _parse_nuclide_symbol(symbol: str) -> Tuple[str, int, int]:
@@ -422,24 +425,28 @@ def get_decay_info_from_paceENSDF(
 
     for k, v in beta_data.items():
         for state in v:
-            # Use ground state branch (level_index == 0) for forbiddenness
+            # Extract per-branch forbiddenness code
-            if state[0] == 0 and len(state) > 15 and state[15] is not None:
+            branch_transition_type = "A"  # default
+            if len(state) > 15 and state[15] is not None:
                 code = str(state[15])
                 if code in FORBIDDENNESS_MAP:
-                    forbidden_code = code
+                    branch_transition_type = FORBIDDENNESS_MAP[code]
-                    transition_type = FORBIDDENNESS_MAP[code]
+
-            # Fallback: use first branch's forbiddenness if no ground state
+            # Use ground state branch (level_index == 0) for overall forbiddenness
-            elif transition_type == "A" and len(state) > 15 and state[15] is not None:
+            if state[0] == 0:
-                code = str(state[15])
+                forbidden_code = (
-                if code in FORBIDDENNESS_MAP:
+                    str(state[15])
-                    forbidden_code = code
+                    if len(state) > 15 and state[15] is not None
-                    transition_type = FORBIDDENNESS_MAP[code]
+                    else "0A"
+                )
+                transition_type = branch_transition_type
 
             branch = BranchInfo(
                 level_index=int(state[0]),
                 level_energy_keV=float(state[1]),
                 intensity=float(state[8]),
                 log_ft=float(state[11]) if state[11] is not None else None,
+                transition_type=branch_transition_type,
             )
             branches.append(branch)
 

beta_spectrum/spectrum.py

@@ -29,6 +29,40 @@ if TYPE_CHECKING:
     from beta_spectrum.components.detector_response import DetectorResponse
 
 
+@dataclass
+class BranchConfig:
+    """
+    Configuration for a single decay branch.
+
+    Each branch corresponds to a decay to a specific daughter nuclear state,
+    with its own endpoint energy (Q_value - excitation_energy) and transition type.
+
+    Attributes
+    ----------
+    endpoint_MeV : float
+        Effective endpoint energy for this branch (Q_value - excitation_energy).
+    transition_type : str
+        Forbiddenness classification for this branch (A, F1, F1U, F2, ...).
+    intensity : float
+        Branch intensity as a fraction of total decays (0.0–1.0).
+        Intensities across all branches should sum to 1.0 (100%).
+    """
+
+    endpoint_MeV: float
+    transition_type: str = "A"
+    intensity: float = 1.0
+
+    def __post_init__(self) -> None:
+        if self.intensity < 0 or self.intensity > 1.0:
+            raise ValueError(
+                f"Branch intensity must be in [0, 1], got {self.intensity}"
+            )
+        if self.endpoint_MeV <= 0:
+            raise ValueError(
+                f"Branch endpoint_MeV must be positive, got {self.endpoint_MeV}"
+            )
+
+
 @dataclass
 class SpectrumConfig:
     """
@@ -37,6 +71,10 @@ class SpectrumConfig:
     Supports declarative detector response specification:
     set detector_model, detector_sigma_a_keV, and related parameters
     to enable automatic detector smearing via convolve_with_detector().
+
+    Supports multi-branch decay via the `branches` field. When branches are
+    provided, the total spectrum is the intensity-weighted sum of individual
+    branch spectra. Each branch has its own endpoint energy and transition type.
     """
 
     Z_parent: int
@@ -69,34 +107,67 @@ class SpectrumConfig:
         0.35,
     )  # in m_e units (total energy)
 
+    # Multi-branch support
+    branches: Optional[List[BranchConfig]] = None  # None → single-branch mode
+    intensity_cutoff: float = (
+        0.0  # Ignore branches below this fraction (default: 0.0 = all)
+    )
+
 
 class BetaSpectrum:
     """
     Precise beta spectrum calculator.
 
-    Combines basic spectrum shape (fermi * phase_space) with mutplicative corrections.
+    Supports both single-branch and multi-branch decay modes.
+
+    In single-branch mode (branches=None), the spectrum is the product of
+    all enabled correction factors.
+
+    In multi-branch mode (branches provided), each branch is calculated
+    independently with its own endpoint energy and transition type, then
+    summed with intensity weighting.
 
     Parameters
     ----------
-    components : List[SpectrumComponent]
+    components : List[SpectrumComponent], optional
-        List of enabled spectral components.
+        List of enabled universal spectral components (Fermi, screening, etc.).
+        Only used in multi-branch mode or for universal components.
+    branch_spectra : List[BetaSpectrum], optional
+        List of per-branch spectrum calculators (each with its own PhaseSpace, etc.).
+        Only used in multi-branch mode.
+    branches : List[BranchConfig], optional
+        Branch configuration metadata (intensities, endpoint energies).
     logger : logging.Logger, optional
         Logger for debug/info output.
     """
 
     def __init__(
         self,
-        components: List[SpectrumComponent],
+        components: Optional[List[SpectrumComponent]] = None,
+        branch_spectra: Optional[List[BetaSpectrum]] = None,
+        branches: Optional[List[BranchConfig]] = None,
         logger: Optional[logging.Logger] = None,
     ):
-        self.components = components
+        self.components = components or []
-        self._component_names = [self._get_component_name(c) for c in components]
+        self._component_names = [self._get_component_name(c) for c in self.components]
+        self.branch_spectra = branch_spectra or []
+        self.branches = branches or []
         self._logger = logger
+        self._is_multi = len(self.branch_spectra) > 0
 
     def __call__(self, W: np.ndarray) -> np.ndarray:
         """
-        Calculate total spectrum as a product of enabled components
+        Calculate spectrum.
+
+        In single-branch mode: product of all components.
+        In multi-branch mode: intensity-weighted sum of branch spectra.
         """
+        if self._is_multi:
+            return self._multi_branch_call(W)
+        return self._single_branch_call(W)
+
+    def _single_branch_call(self, W: np.ndarray) -> np.ndarray:
+        """Calculate spectrum as product of components (single branch)."""
         result: np.ndarray = np.ones_like(W, dtype=float)
         for comp in self.components:
             if self._logger:
@@ -106,21 +177,74 @@ class BetaSpectrum:
             result = result * comp(W)
         return result.astype(np.float64)
 
+    def _multi_branch_call(self, W: np.ndarray) -> np.ndarray:
+        """Calculate total spectrum as intensity-weighted sum of branches.
+
+        The total is: (universal_components_product) * sum(intensity_i * branch_spectrum_i)
+        """
+        # Sum of intensity-weighted branch spectra
+        branch_sum = np.zeros_like(W, dtype=float)
+        for i, branch_spec in enumerate(self.branch_spectra):
+            if self._logger:
+                self._logger.debug(
+                    "Evaluating branch %d/%d (endpoint=%.3f MeV, intensity=%.4f)",
+                    i + 1,
+                    len(self.branch_spectra),
+                    self.branches[i].endpoint_MeV,
+                    self.branches[i].intensity,
+                )
+            branch_spectrum = branch_spec(W)
+            # Mask out energies beyond this branch's endpoint
+            W0_branch = T_to_W(self.branches[i].endpoint_MeV)
+            mask = W <= W0_branch
+            branch_sum += self.branches[i].intensity * np.where(
+                mask, branch_spectrum, 0.0
+            )
+
+        # Multiply by universal components
+        result: np.ndarray = np.ones_like(W, dtype=float)
+        for comp in self.components:
+            result = result * comp(W)
+
+        return np.asarray((result * branch_sum).astype(np.float64), dtype=np.float64)
+
     def _get_component_name(self, comp: SpectrumComponent) -> str:
-        """
+        """Get human-readable component name."""
-        Get human-readable component names
-        """
         class_name = comp.__class__.__name__
         return class_name.replace("Correction", "").replace("Function", "")
 
     def calculate_components(self, W: np.ndarray) -> Dict[str, np.ndarray]:
         """
-        Calculate individual components for debugging and analysis
+        Calculate individual components.
+
+        In single-branch mode: returns per-component values.
+        In multi-branch mode: returns universal components + per-branch
+        branch-dependent components (PhaseSpace, Radiative).
         """
-        components_dict = {}
+        if self._is_multi:
+            return self._multi_branch_components(W)
+
+        result: Dict[str, np.ndarray] = {}
         for name, comp in zip(self._component_names, self.components):
-            components_dict[name] = comp(W)
+            result[name] = comp(W)
-        return components_dict
+        return result
+
+    def _multi_branch_components(self, W: np.ndarray) -> Dict[str, np.ndarray]:
+        """Return universal + per-branch components for multi-branch mode."""
+        result: Dict[str, np.ndarray] = {}
+
+        # Universal components (same for all branches)
+        for name, comp in zip(self._component_names, self.components):
+            result[name] = comp(W)
+
+        # Per-branch components (branch-dependent)
+        for i, branch_spec in enumerate(self.branch_spectra):
+            branch_name = f"branch_{i}"
+            branch_components = branch_spec.calculate_components(W)
+            for comp_name, values in branch_components.items():
+                result[f"{branch_name}.{comp_name}"] = values
+
+        return result
 
     @classmethod
     def from_config(
@@ -129,6 +253,9 @@ class BetaSpectrum:
         """
         Create a BetaSpectrum from configuration.
 
+        In single-branch mode: creates a single spectrum calculator.
+        In multi-branch mode: creates per-branch calculators + universal components.
+
         Parameters
         ----------
         config : SpectrumConfig
@@ -136,6 +263,15 @@ class BetaSpectrum:
         logger : logging.Logger, optional
             Logger for progress output.
         """
+        if config.branches is not None and len(config.branches) > 1:
+            return cls._from_config_multi(config, logger)
+        return cls._from_config_single(config, logger)
+
+    @classmethod
+    def _from_config_single(
+        cls, config: SpectrumConfig, logger: Optional[logging.Logger] = None
+    ) -> BetaSpectrum:
+        """Create single-branch spectrum calculator."""
         components: List[SpectrumComponent] = []
         W0 = float(T_to_W(config.endpoint_MeV))
 
@@ -185,17 +321,163 @@ class BetaSpectrum:
 
         return cls(components, logger=logger)
 
+    @classmethod
+    def _from_config_multi(
+        cls, config: SpectrumConfig, logger: Optional[logging.Logger] = None
+    ) -> BetaSpectrum:
+        """Create multi-branch spectrum calculator.
+
+        Creates per-branch calculators (each with its own PhaseSpace, Radiative, etc.)
+        and universal components (Fermi, screening, finite_size, etc.) shared across
+        all branches.
+
+        Note: Radiative correction is branch-dependent (depends on W0), so it goes
+        into each branch's calculator, not the universal components.
+        """
+        if logger:
+            logger.info("Multi-branch mode: %d branches", len(config.branches))
+
+        # Create per-branch calculators
+        branch_spectra: List[BetaSpectrum] = []
+        for i, branch in enumerate(config.branches):
+            if logger:
+                logger.info(
+                    "  Branch %d: endpoint=%.3f MeV, transition=%s, intensity=%.4f",
+                    i + 1,
+                    branch.endpoint_MeV,
+                    branch.transition_type,
+                    branch.intensity,
+                )
+
+            # Create a per-branch config with the branch's endpoint and transition type
+            branch_config = SpectrumConfig(
+                Z_parent=config.Z_parent,
+                Z_daughter=config.Z_daughter,
+                A_number=config.A_number,
+                endpoint_MeV=branch.endpoint_MeV,
+                transition_type=branch.transition_type,
+                e_step_MeV=config.e_step_MeV,
+                # Branch-dependent components
+                use_phase_space=config.use_phase_space,
+                use_radiative=config.use_radiative,
+                # Universal components disabled (handled separately)
+                use_fermi=False,
+                use_screening=False,
+                use_finite_size=False,
+                use_charge_dist=False,
+                use_exchange=False,
+            )
+
+            branch_spectrum = cls._from_config_single(branch_config, logger)
+            branch_spectra.append(branch_spectrum)
+
+        # Create universal components (shared across all branches)
+        # Note: Radiative is NOT included here because it depends on W0
+        universal_components: List[SpectrumComponent] = []
+
+        if config.use_fermi:
+            universal_components.append(
+                FermiFunction(Z=config.Z_daughter, A=config.A_number, logger=logger)
+            )
+
+        if config.use_finite_size:
+            universal_components.append(
+                FiniteSizeL0(Z=config.Z_daughter, A=config.A_number, logger=logger)
+            )
+
+        if config.use_charge_dist:
+            universal_components.append(
+                ChargeDistributionU(
+                    Z=config.Z_daughter, A=config.A_number, logger=logger
+                )
+            )
+
+        if config.use_screening:
+            universal_components.append(
+                ScreeningCorrection(
+                    FermiFunction(Z=config.Z_parent, A=config.A_number, logger=logger),
+                    logger=logger,
+                )
+            )
+
+        if config.use_exchange:
+            universal_components.append(
+                ExchangeCorrection(Z=config.Z_parent, logger=logger)
+            )
+
+        return cls(
+            components=universal_components,
+            branch_spectra=branch_spectra,
+            branches=list(config.branches),
+            logger=logger,
+        )
+
     def get_energy_grid(self, config: SpectrumConfig) -> Tuple[np.ndarray, np.ndarray]:
         """
         Generate array of energy values for spectrum calculation.
+
+        In multi-branch mode, uses the maximum endpoint across all branches.
         """
-        kinetic_MeV = np.arange(
+        if self._is_multi and self.branches:
-            config.e_step_MeV, config.endpoint_MeV, config.e_step_MeV
+            max_endpoint = max(b.endpoint_MeV for b in self.branches)
-        )
+        else:
+            max_endpoint = config.endpoint_MeV
+
+        kinetic_MeV = np.arange(config.e_step_MeV, max_endpoint, config.e_step_MeV)
         W = T_to_W(kinetic_MeV)
 
         return W, kinetic_MeV.astype(np.float64)
 
+    def get_branch_spectra(self, W: np.ndarray) -> List[np.ndarray]:
+        """
+        Get individual branch spectra for multi-branch mode.
+
+        Each branch spectrum is masked to zero beyond its endpoint.
+
+        Parameters
+        ----------
+        W : np.ndarray
+            Energy grid in m_e units.
+
+        Returns
+        -------
+        List[np.ndarray]
+            Per-branch spectrum arrays (unweighted). Empty in single-branch mode.
+        """
+        if not self._is_multi:
+            return []
+        result = []
+        for i, branch_spec in enumerate(self.branch_spectra):
+            branch_spectrum = branch_spec(W)
+            # Use np.maximum to smoothly approach zero at endpoint,
+            # and clamp negative values beyond endpoint to zero
+            result.append(np.maximum(0.0, branch_spectrum))
+        return result
+
+    def get_branch_normalized_spectra(self, W: np.ndarray) -> List[np.ndarray]:
+        """
+        Get intensity-weighted branch spectra for multi-branch mode.
+
+        Each branch spectrum is masked to zero beyond its endpoint.
+
+        Parameters
+        ----------
+        W : np.ndarray
+            Energy grid in m_e units.
+
+        Returns
+        -------
+        List[np.ndarray]
+            Per-branch spectrum arrays weighted by intensity. Empty in single-branch mode.
+        """
+        if not self._is_multi:
+            return []
+        result = []
+        for i, branch_spec in enumerate(self.branch_spectra):
+            branch_spectrum = branch_spec(W)
+            result.append(self.branches[i].intensity * np.maximum(0.0, branch_spectrum))
+        return result
+
     @staticmethod
     def create_detector_from_config(
         config: SpectrumConfig,
@@ -353,7 +635,9 @@ class BetaSpectrumAnalyzer:
         self.config = config
         self.W, self.energies_MeV = spectrum.get_energy_grid(config)
         self._components_cache: Optional[Dict[str, np.ndarray]] = None
+        self._branch_spectra_cache: Optional[List[np.ndarray]] = None
         self._logger = logger
+        self._is_multi = spectrum._is_multi
 
     @property
     def components(self) -> Dict[str, np.ndarray]:
@@ -362,6 +646,16 @@ class BetaSpectrumAnalyzer:
         assert self._components_cache is not None
         return self._components_cache
 
+    @property
+    def branch_spectra(self) -> List[np.ndarray]:
+        """Intensity-weighted branch spectra."""
+        if self._branch_spectra_cache is None:
+            self._branch_spectra_cache = self.spectrum.get_branch_normalized_spectra(
+                self.W
+            )
+        assert self._branch_spectra_cache is not None
+        return self._branch_spectra_cache
+
     def total_spectrum(self, normalize: bool = True) -> np.ndarray:
         """
         Calculate the total_spectrum
@@ -425,23 +719,39 @@ class BetaSpectrumAnalyzer:
         """
         Create visualization of the spectrum and all correction factors.
 
+        In single-branch mode: standard 4-panel debug view or single spectrum plot.
+        In multi-branch mode: spectrum plot with individual branches shown,
+        and debug view with vertical layout for branch comparison.
+
         Parameters
         ----------
         save_path : str, optional
             Path to save the figure.
         show_components : bool
-            If True, show full debug view (4-panel with all corrections).
+            If True, show full debug view.
             If False, show only the spectrum plot with nuclear data header.
+            If True, use logarithmic y-scale. Set to False for
+            linear scale, which is better when branch intensities are similar.
         """
         total = self.total_spectrum(normalize=True)
         components = self.components
         commit = get_git_short_hash(6)
         timestamp = datetime.now(timezone.utc).strftime("%Y-%m-%d %H:%M:%S UTC")
 
-        if show_components:
+        if self._is_multi:
-            self._plot_debug_view(total, components, save_path, commit, timestamp)
+            if show_components:
+                self._plot_multi_branch_debug(
+                    total, components, save_path, commit, timestamp
+                )
+            else:
+                self._plot_multi_branch_spectra(
+                    total, components, save_path, commit, timestamp
+                )
         else:
-            self._plot_spectrum_only(total, save_path, commit, timestamp)
+            if show_components:
+                self._plot_debug_view(total, components, save_path, commit, timestamp)
+            else:
+                self._plot_spectrum_only(total, save_path, commit, timestamp)
 
     def _add_id_textbox(self, ax: "Axes", commit: str, timestamp: str) -> None:
         """Add ID text box with commit hash and timestamp to plot."""
@@ -780,6 +1090,302 @@ class BetaSpectrumAnalyzer:
         ]
         return symbols[Z] if Z < len(symbols) else f"Z{Z}"
 
+    # ---------------------------------------------------------------------------
+    # Multi-branch plotting
+    # ---------------------------------------------------------------------------
+
+    def _get_branch_label(self, i: int) -> str:
+        """Generate branch label for legend."""
+        if i >= len(self.spectrum.branches):
+            return f"Branch {i+1}"
+        branch = self.spectrum.branches[i]
+        return f"Branch {i+1}: {branch.transition_type}, E₀={branch.endpoint_MeV*1000:.1f} keV"
+
+    def _plot_multi_branch_spectra(
+        self,
+        total: np.ndarray,
+        components: Dict[str, np.ndarray],
+        save_path: Optional[str],
+        commit: str,
+        timestamp: str,
+    ) -> None:
+        """Plot total spectrum with individual branch decomposition.
+
+        All curves use raw (un-normalized) values so that the total is
+        the arithmetic sum of branch contributions, and branch amplitudes
+        are proportional to their intensities.
+        """
+        fig, ax = plt.subplots(figsize=(12, 7))
+
+        branch_colors = [
+            "#1f77b4",
+            "#ff7f0e",
+            "#2ca02c",
+            "#d62728",
+            "#9467bd",
+            "#8c564b",
+            "#e377c2",
+            "#7f7f7f",
+            "#bcbd22",
+            "#17becf",
+        ]
+
+        # Compute branch contributions (universal * intensity * branch_spectrum)
+        branch_contribs = []
+        for i, branch_spec in enumerate(self.spectrum.branch_spectra):
+            universal_product = np.ones_like(self.W)
+            for name in [
+                "Fermi",
+                "Screening",
+                "FiniteSizeL0",
+                "ChargeDistributionU",
+                "Exchange",
+            ]:
+                if name in components:
+                    universal_product *= components[name]
+            bs = self.spectrum.get_branch_spectra(self.W)
+            contrib = universal_product * self.spectrum.branches[i].intensity * bs[i]
+            branch_contribs.append(contrib)
+
+        # Total = arithmetic sum of all branches (raw, un-normalized)
+        total_raw = sum(branch_contribs)
+
+        ax.plot(
+            self.energies_MeV,
+            total_raw,
+            "k-",
+            lw=2.5,
+            label="Total spectrum",
+            zorder=10,
+        )
+
+        # Plot each branch contribution (raw, un-normalized)
+        for i, (contrib, branch) in enumerate(
+            zip(branch_contribs, self.spectrum.branches)
+        ):
+            color = branch_colors[i % len(branch_colors)]
+            ax.plot(
+                self.energies_MeV,
+                contrib,
+                color=color,
+                lw=1.5,
+                alpha=0.7,
+                label=self._get_branch_label(i),
+            )
+
+        ax.set_xlabel(r"Electron kinetic energy $E$ [MeV]", fontsize=12)
+        ax.set_ylabel("Normalized Counts", fontsize=12)
+        parent = self._element_symbol(self.config.Z_parent)
+        daughter = self._element_symbol(self.config.Z_daughter)
+        ax.set_title(
+            f"Beta-decay: {parent}{self.config.A_number} -> {daughter}{self.config.A_number} — {len(self.spectrum.branches)} branches",
+            fontsize=14,
+        )
+        ax.set_yscale("log")
+        ax.grid(True, alpha=0.3)
+        ax.legend(loc="best", fontsize=10, ncol=min(2, len(self.spectrum.branches) + 1))
+
+        self._add_nuclear_data_header(ax)
+        self._add_id_textbox(ax, commit, timestamp)
+
+        plt.tight_layout()
+
+        if save_path:
+            fig.savefig(save_path, dpi=150, bbox_inches="tight")
+            print(f"Figure saved to {save_path}")
+
+    def _plot_multi_branch_debug(
+        self,
+        total: np.ndarray,
+        components: Dict[str, np.ndarray],
+        save_path: Optional[str],
+        commit: str,
+        timestamp: str,
+    ) -> None:
+        """Debug view for multi-branch: vertical layout with 4 panels.
+
+        Panel 1: Total spectrum + intensity-scaled branch spectra (log scale)
+        Panel 2: Fermi function + all phase spaces (shape comparison, no intensity)
+        Panel 3: Universal W0-independent corrections (Fermi, Screening, etc.)
+        Panel 4: W0-dependent corrections (Radiative for each branch)
+        """
+        n_branches = len(self.spectrum.branches)
+
+        fig, axes = plt.subplots(4, 1, figsize=(12, 16), sharex=True)
+
+        branch_colors = [
+            "#1f77b4",
+            "#ff7f0e",
+            "#2ca02c",
+            "#d62728",
+            "#9467bd",
+            "#8c564b",
+            "#e377c2",
+            "#7f7f7f",
+            "#bcbd22",
+            "#17becf",
+        ]
+        parent = self._element_symbol(self.config.Z_parent)
+        daughter = self._element_symbol(self.config.Z_daughter)
+
+        # Panel 1: Total spectrum + per-branch spectra
+        # All curves use raw (un-normalized) values so total = sum of branches
+        branch_contribs = []
+        for i, branch_spec in enumerate(self.spectrum.branch_spectra):
+            universal_product = np.ones_like(self.W)
+            for name in [
+                "Fermi",
+                "Screening",
+                "FiniteSizeL0",
+                "ChargeDistributionU",
+                "Exchange",
+            ]:
+                if name in components:
+                    universal_product *= components[name]
+            bs = self.spectrum.get_branch_spectra(self.W)
+            contrib = universal_product * self.spectrum.branches[i].intensity * bs[i]
+            branch_contribs.append(contrib)
+
+        total_raw = sum(branch_contribs)
+
+        ax = axes[0]
+        ax.plot(
+            self.energies_MeV,
+            total_raw,
+            "k-",
+            lw=2,
+            label="Total spectrum",
+            zorder=10,
+        )
+
+        for i, (contrib, branch) in enumerate(
+            zip(branch_contribs, self.spectrum.branches)
+        ):
+            color = branch_colors[i % len(branch_colors)]
+            ax.plot(
+                self.energies_MeV,
+                contrib,
+                color=color,
+                lw=1.5,
+                alpha=0.7,
+                label=f"Branch {i+1}: {branch.transition_type}, "
+                f"E₀={branch.endpoint_MeV*1000:.1f} keV",
+            )
+        ax.set_yscale("log")
+        ax.set_ylabel("Normalized Counts", fontsize=10)
+        ax.set_title(
+            f"Multi-branch: {parent}{self.config.A_number} -> {daughter}{self.config.A_number} ({n_branches} branches)",
+            fontsize=12,
+            fontweight="bold",
+        )
+        ax.grid(True, alpha=0.3)
+        ax.legend(fontsize=8, loc="best", ncol=min(2, n_branches + 1))
+
+        # Panel 2: Fermi function + all phase spaces (shape comparison)
+        # Phase space is masked beyond endpoint before normalizing to avoid
+        # negative values from the sqrt term
+        ax = axes[1]
+        if "Fermi" in components:
+            fermi_vals = components["Fermi"]
+            fermi_scaled = fermi_vals / np.max(fermi_vals)
+            ax.plot(
+                self.energies_MeV,
+                fermi_scaled,
+                "r-",
+                lw=1.5,
+                label="Fermi (normalized)",
+                alpha=0.8,
+            )
+        for i, (branch_spec, branch) in enumerate(
+            zip(self.spectrum.branch_spectra, self.spectrum.branches)
+        ):
+            color = branch_colors[i % len(branch_colors)]
+            if f"branch_{i}.PhaseSpace" in components:
+                ps_vals = components[f"branch_{i}.PhaseSpace"]
+                # Mask beyond endpoint to avoid negative values
+                W0_branch = T_to_W(branch.endpoint_MeV)
+                mask = self.W <= W0_branch
+                ps_masked = np.where(mask, ps_vals, 0.0)
+                ps_max = np.max(ps_masked)
+                if ps_max > 0:
+                    ps_scaled = ps_masked / ps_max
+                else:
+                    ps_scaled = ps_masked
+                ax.plot(
+                    self.energies_MeV,
+                    ps_scaled,
+                    color=color,
+                    lw=1.5,
+                    alpha=0.7,
+                    label=f"PS Branch {i+1} (normalized)",
+                )
+        ax.set_ylabel("Normalized Factor", fontsize=10)
+        ax.set_title(
+            "Fermi Function & Phase Space Shapes", fontsize=11, fontweight="bold"
+        )
+        ax.grid(True, alpha=0.3)
+        ax.legend(fontsize=8, loc="best")
+
+        # Panel 3: Universal W0-independent corrections
+        ax = axes[2]
+        universal_components = [
+            "Screening",
+            "FiniteSizeL0",
+            "ChargeDistributionU",
+            "Exchange",
+        ]
+        colors = ["green", "orange", "purple", "brown"]
+        for name, color in zip(universal_components, colors):
+            if name in components:
+                values = components[name]
+                values_scaled = values / np.max(values)
+                ax.plot(
+                    self.energies_MeV,
+                    values_scaled,
+                    color=color,
+                    lw=1.5,
+                    label=f"{name} (normalized)",
+                    alpha=0.8,
+                )
+        ax.set_ylabel("Normalized Factor", fontsize=10)
+        ax.set_title(
+            "Universal W0-Independent Corrections", fontsize=11, fontweight="bold"
+        )
+        ax.grid(True, alpha=0.3)
+        ax.legend(fontsize=8, loc="best")
+
+        # Panel 4: W0-dependent corrections (Radiative for each branch)
+        ax = axes[3]
+        for i, (branch_spec, branch) in enumerate(
+            zip(self.spectrum.branch_spectra, self.spectrum.branches)
+        ):
+            color = branch_colors[i % len(branch_colors)]
+            if f"branch_{i}.Radiative" in components:
+                rad_vals = components[f"branch_{i}.Radiative"]
+                ax.plot(
+                    self.energies_MeV,
+                    rad_vals,
+                    color=color,
+                    lw=1.5,
+                    alpha=0.7,
+                    label=f"Branch {i+1} Radiative (E₀={branch.endpoint_MeV*1000:.1f} keV)",
+                )
+        ax.set_ylabel("Radiative Correction", fontsize=10)
+        ax.set_xlabel(r"Electron kinetic energy $E$ [MeV]", fontsize=10)
+        ax.set_title(
+            "W0-Dependent Corrections (Radiative)", fontsize=11, fontweight="bold"
+        )
+        ax.grid(True, alpha=0.3)
+        ax.legend(fontsize=8, loc="best", ncol=min(2, n_branches + 1))
+
+        self._add_id_textbox(axes[3], commit, timestamp)
+
+        plt.tight_layout()
+
+        if save_path:
+            fig.savefig(save_path, dpi=150, bbox_inches="tight")
+            print(f"Figure saved to {save_path}")
+
     def export_to_csv(self, filename: str, source_type: str = "unknown") -> None:
         """
         Export spectrum data to CSV with metadata header.
@@ -934,6 +1540,17 @@ class BetaSpectrumAnalyzer:
             f"# e_step: {self.config.e_step_MeV:.4f} MeV",
         ]
 
+        if self._is_multi and self.spectrum.branches:
+            header_lines.append(f"# branches: {len(self.spectrum.branches)}")
+            for i, branch in enumerate(self.spectrum.branches):
+                header_lines.append(
+                    f"#   branch_{i+1}: E₀={branch.endpoint_MeV*1000:.1f} keV, "
+                    f"transition={branch.transition_type}, "
+                    f"intensity={branch.intensity:.4f}"
+                )
+        else:
+            header_lines.append("# branches: 1 (single)")
+
         if self.config.use_detector_response:
             header_lines.append(
                 f"# detector: {self.config.detector_model} "
@@ -952,8 +1569,37 @@ class BetaSpectrumAnalyzer:
             "energy_MeV": self.energies_MeV,
             "spectrum": total,
         }
-        for name, values in components.items():
+
-            data[name] = values
+        if self._is_multi and self.spectrum.branches:
+            # Multi-branch mode: total spectrum + per-branch spectra + components
+            for i, (branch, branch_spec) in enumerate(
+                zip(self.spectrum.branches, self.spectrum.branch_spectra)
+            ):
+                data[f"branch_{i+1}_spectrum"] = branch_spec(self.W)
+                data[f"branch_{i+1}_intensity"] = np.full_like(
+                    self.energies_MeV, branch.intensity
+                )
+                data[f"branch_{i+1}_endpoint"] = np.full_like(
+                    self.energies_MeV, branch.endpoint_MeV
+                )
+                data[f"branch_{i+1}_transition"] = np.full_like(
+                    self.energies_MeV, branch.transition_type, dtype=object
+                )
+
+            # Per-branch components (PhaseSpace, Radiative, etc.)
+            for i, branch_spec in enumerate(self.spectrum.branch_spectra):
+                branch_comps = branch_spec.calculate_components(self.W)
+                for comp_name, values in branch_comps.items():
+                    data[f"branch_{i+1}_{comp_name}"] = values
+
+            # Universal components only (skip per-branch keys)
+            for name, values in components.items():
+                if not name.startswith("branch_"):
+                    data[f"universal_{name}"] = values
+        else:
+            # Single-branch mode: standard export
+            for name, values in components.items():
+                data[name] = values
 
         df = pd.DataFrame(data)
         df.to_csv(filename, index=False, float_format="%.4e")
@@ -965,7 +1611,8 @@ class BetaSpectrumAnalyzer:
             f.write("\n".join(header_lines) + "\n" + content)
 
         if self._logger:
-            self._logger.info("CSV exported to %s with metadata header", filename)
+            mode = "multi-branch" if self._is_multi else "single-branch"
+            self._logger.info("CSV exported to %s (%s mode)", filename, mode)
 
     def get_data(self) -> Dict[str, Any]:
         """Get all numerical data for custom analysis."""

tests/test_multi_branch.py

@@ -0,0 +1,775 @@
+"""Tests for multi-branch decay support.
+
+Why test this?
+--------------
+Multi-branch decay is a critical feature for realistic nuclear decay analysis.
+Many isotopes decay through multiple channels to different daughter states,
+each with its own endpoint energy, transition type, and intensity.
+
+The multi-branch mode is the DEFAULT — even single-branch decays are handled
+through the same code path. The only difference in behavior is in plot generation
+(where a single branch doesn't need separate branch plots).
+
+We test:
+
+ 1. **BranchConfig validation**: intensity bounds, positive endpoint.
+ 2. **Single-branch via multi-branch path**: one branch should give identical
+    results to the old single-branch behavior.
+ 3. **Multi-branch spectrum calculation**: weighted sum produces valid output.
+ 4. **Intensity weighting**: sum of intensities = 1.0, weighted sum = total.
+ 5. **Per-branch spectra**: each branch has correct endpoint and transition type.
+ 6. **Component extraction**: universal + per-branch components correctly separated.
+ 7. **CSV export**: multi-branch CSV has correct column structure.
+ 8. **Energy grid**: uses max endpoint across all branches.
+ 9. **Detector response**: multi-branch spectra are convolvable.
+ 10. **BranchConfig dataclass**: basic construction and validation.
+"""
+
+import numpy as np
+import pytest
+
+from beta_spectrum import (
+    BetaSpectrum,
+    BetaSpectrumAnalyzer,
+    BranchConfig,
+    SpectrumConfig,
+)
+from beta_spectrum.utils import T_to_W
+
+# ---------------------------------------------------------------------------
+# BranchConfig validation
+# ---------------------------------------------------------------------------
+
+
+class TestBranchConfigValidation:
+    """Test BranchConfig dataclass validation."""
+
+    def test_valid_branch(self):
+        """Valid branch config should construct without error."""
+        branch = BranchConfig(
+            endpoint_MeV=0.2,
+            transition_type="A",
+            intensity=0.5,
+        )
+        assert branch.endpoint_MeV == 0.2
+        assert branch.transition_type == "A"
+        assert branch.intensity == 0.5
+
+    def test_intensity_out_of_range_low(self):
+        """Intensity below 0 should raise ValueError."""
+        with pytest.raises(ValueError, match="intensity must be in"):
+            BranchConfig(endpoint_MeV=0.2, intensity=-0.1)
+
+    def test_intensity_out_of_range_high(self):
+        """Intensity above 1.0 should raise ValueError."""
+        with pytest.raises(ValueError, match="intensity must be in"):
+            BranchConfig(endpoint_MeV=0.2, intensity=1.1)
+
+    def test_zero_endpoint(self):
+        """Zero endpoint should raise ValueError."""
+        with pytest.raises(ValueError, match="endpoint_MeV must be positive"):
+            BranchConfig(endpoint_MeV=0.0, intensity=0.5)
+
+    def test_negative_endpoint(self):
+        """Negative endpoint should raise ValueError."""
+        with pytest.raises(ValueError, match="endpoint_MeV must be positive"):
+            BranchConfig(endpoint_MeV=-0.1, intensity=0.5)
+
+    def test_default_transition_type(self):
+        """Transition type defaults to 'A'."""
+        branch = BranchConfig(endpoint_MeV=0.2, intensity=0.5)
+        assert branch.transition_type == "A"
+
+
+# ---------------------------------------------------------------------------
+# Single-branch via multi-branch path (backward compatibility)
+# ---------------------------------------------------------------------------
+
+
+class TestSingleBranchViaMultiBranch:
+    """Verify single-branch mode works through the multi-branch path."""
+
+    def test_one_branch_produces_valid_spectrum(self):
+        """One branch should produce valid spectrum identical to old behavior."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.294,
+            branches=[
+                BranchConfig(
+                    endpoint_MeV=0.294,
+                    transition_type="F2",
+                    intensity=1.0,
+                )
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, E = spectrum.get_energy_grid(config)
+        values = spectrum(W)
+
+        assert len(values) == len(E)
+        assert np.all(np.isfinite(values))
+        assert np.all(values >= 0), "Spectrum should be non-negative"
+
+    def test_one_branch_has_no_branch_spectra(self):
+        """Single branch should not produce separate branch plots."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.294,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=1.0),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, _ = spectrum.get_energy_grid(config)
+
+        branch_spectra = spectrum.get_branch_spectra(W)
+        # Single branch: no separate branch spectra needed
+        assert len(branch_spectra) == 0
+
+    def test_one_branch_analyzer_uses_normal_plot(self):
+        """Single branch analyzer should not have multi-branch mode."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.294,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=1.0),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+
+        assert not analyzer._is_multi
+
+
+# ---------------------------------------------------------------------------
+# Multi-branch spectrum calculation
+# ---------------------------------------------------------------------------
+
+
+class TestMultiBranchSpectrum:
+    """Test multi-branch spectrum calculation."""
+
+    def test_two_branches_basic(self):
+        """Two branches should produce valid weighted sum."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+
+        assert spectrum._is_multi
+        assert len(spectrum.branch_spectra) == 2
+        assert len(spectrum.branches) == 2
+
+        W, E = spectrum.get_energy_grid(config)
+        values = spectrum(W)
+
+        assert len(values) == len(E)
+        assert np.all(np.isfinite(values))
+        assert np.all(values >= 0), "Spectrum should be non-negative"
+
+    def test_three_branches(self):
+        """Three branches should produce valid weighted sum."""
+        config = SpectrumConfig(
+            Z_parent=27,
+            Z_daughter=28,
+            A_number=60,
+            endpoint_MeV=0.320,
+            branches=[
+                BranchConfig(endpoint_MeV=0.310, transition_type="A", intensity=0.5),
+                BranchConfig(endpoint_MeV=0.200, transition_type="F1", intensity=0.3),
+                BranchConfig(endpoint_MeV=0.100, transition_type="A", intensity=0.2),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+
+        W, E = spectrum.get_energy_grid(config)
+        values = spectrum(W)
+
+        assert spectrum._is_multi
+        assert len(values) == len(E)
+        assert np.all(np.isfinite(values))
+        assert np.all(values >= 0)
+
+    def test_energy_grid_uses_max_endpoint(self):
+        """Energy grid should extend to the maximum branch endpoint."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            e_step_MeV=0.001,
+            branches=[
+                BranchConfig(endpoint_MeV=0.300, transition_type="A", intensity=0.5),
+                BranchConfig(endpoint_MeV=0.200, transition_type="F2", intensity=0.5),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, E = spectrum.get_energy_grid(config)
+
+        # Grid should extend to max endpoint (0.300 MeV)
+        assert E[-1] >= 0.299, f"Grid should reach ~0.300 MeV, got {E[-1]}"
+
+    def test_branch_spectra_property(self):
+        """get_branch_spectra() should return per-branch unweighted spectra."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, _ = spectrum.get_energy_grid(config)
+
+        branch_spectra = spectrum.get_branch_spectra(W)
+        assert len(branch_spectra) == 2
+        assert len(branch_spectra[0]) == len(W)
+        assert len(branch_spectra[1]) == len(W)
+
+    def test_branch_normalized_spectra(self):
+        """get_branch_normalized_spectra() should return intensity-weighted spectra."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            # Disable universal components to simplify
+            use_fermi=False,
+            use_screening=False,
+            use_finite_size=False,
+            use_charge_dist=False,
+            use_exchange=False,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, _ = spectrum.get_energy_grid(config)
+
+        weighted = spectrum.get_branch_normalized_spectra(W)
+        assert len(weighted) == 2
+
+        # With no universal components, total = sum of weighted branch spectra
+        total_from_branches = weighted[0] + weighted[1]
+        total = spectrum(W)
+        assert np.allclose(total, total_from_branches, rtol=1e-10)
+
+    def test_intensity_normalization(self):
+        """Branch intensities should be used as weights correctly."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            use_phase_space=True,
+            use_fermi=False,
+            use_screening=False,
+            use_finite_size=False,
+            use_charge_dist=False,
+            use_exchange=False,
+            use_radiative=False,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="A", intensity=0.8),
+                BranchConfig(endpoint_MeV=0.294, transition_type="A", intensity=0.2),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, _ = spectrum.get_energy_grid(config)
+
+        # Both branches have same endpoint and transition type,
+        # so phase space should be identical. Total = 0.8*PS + 0.2*PS = PS
+        total = spectrum(W)
+        branch0 = spectrum.branch_spectra[0](W)
+        branch1 = spectrum.branch_spectra[1](W)
+
+        expected = 0.8 * branch0 + 0.2 * branch1
+        assert np.allclose(total, expected, rtol=1e-10)
+
+
+# ---------------------------------------------------------------------------
+# Component extraction in multi-branch mode
+# ---------------------------------------------------------------------------
+
+
+class TestMultiBranchComponents:
+    """Test component extraction for multi-branch mode."""
+
+    def test_universal_components_present(self):
+        """Universal components (Fermi, etc.) should be in component dict."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, _ = spectrum.get_energy_grid(config)
+        components = spectrum.calculate_components(W)
+
+        assert "Fermi" in components, "Universal Fermi component should be present"
+        assert "Screening" in components
+        assert "FiniteSizeL0" in components
+
+    def test_per_branch_components_present(self):
+        """Per-branch components (PhaseSpace, Radiative) should be in component dict."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, _ = spectrum.get_energy_grid(config)
+        components = spectrum.calculate_components(W)
+
+        assert "branch_0.PhaseSpace" in components
+        assert "branch_1.PhaseSpace" in components
+        assert "branch_0.Radiative" in components
+        assert "branch_1.Radiative" in components
+
+    def test_analyzer_components_multi_branch(self):
+        """BetaSpectrumAnalyzer should handle multi-branch components."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+
+        components = analyzer.components
+        assert "Fermi" in components
+        assert "branch_0.PhaseSpace" in components
+
+    def test_analyzer_branch_spectra_property(self):
+        """Analyzer should expose branch spectra property."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+
+        branch_spectra = analyzer.branch_spectra
+        assert len(branch_spectra) == 2
+
+    def test_total_equals_universal_times_branch_sum(self):
+        """Total = universal_components * sum(intensity * branch_spectra)."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, _ = spectrum.get_energy_grid(config)
+
+        total = spectrum(W)
+
+        # Compute expected: universal_components * sum(intensity * branch_spectra)
+        # Each branch spectrum must be masked to zero beyond its endpoint
+        universal_product = np.ones_like(W)
+        for comp in spectrum.components:
+            universal_product *= comp(W)
+
+        branch_sum = np.zeros_like(W)
+        for i, spec in enumerate(spectrum.branch_spectra):
+            W0_branch = T_to_W(spectrum.branches[i].endpoint_MeV)
+            mask = W <= W0_branch
+            branch_sum += spectrum.branches[i].intensity * np.where(mask, spec(W), 0.0)
+
+        expected = universal_product * branch_sum
+        assert np.allclose(total, expected, rtol=1e-10)
+
+
+# ---------------------------------------------------------------------------
+# CSV export in multi-branch mode
+# ---------------------------------------------------------------------------
+
+
+class TestMultiBranchCSVExport:
+    """Test CSV export structure for multi-branch mode."""
+
+    def _write_json(self, data):
+        import json
+        import tempfile
+        from pathlib import Path
+
+        tmp = tempfile.NamedTemporaryFile(mode="w", suffix=".json", delete=False)
+        json.dump(data, tmp)
+        tmp.close()
+        return Path(tmp.name)
+
+    def test_multi_branch_csv_columns(self, tmp_path):
+        """Multi-branch CSV should have total + per-branch + universal columns."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+
+        csv_path = str(tmp_path / "multi_branch.csv")
+        analyzer.export_to_csv(csv_path)
+
+        import pandas as pd
+
+        # Read CSV (skip comment headers)
+        df = pd.read_csv(csv_path, comment="#")
+
+        # Should have total spectrum column
+        assert "spectrum" in df.columns
+
+        # Should have per-branch columns
+        assert "branch_1_spectrum" in df.columns
+        assert "branch_2_spectrum" in df.columns
+
+        # Should have universal columns
+        assert "universal_Fermi" in df.columns
+
+    def test_single_branch_csv_columns(self, tmp_path):
+        """Single-branch CSV should have standard columns."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.294,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=1.0),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        analyzer = BetaSpectrumAnalyzer(spectrum, config)
+
+        csv_path = str(tmp_path / "single_branch.csv")
+        analyzer.export_to_csv(csv_path)
+
+        import pandas as pd
+
+        df = pd.read_csv(csv_path, comment="#")
+
+        assert "spectrum" in df.columns
+        assert "Fermi" in df.columns
+
+
+# ---------------------------------------------------------------------------
+# Multi-branch with different transition types
+# ---------------------------------------------------------------------------
+
+
+class TestMultiBranchTransitionTypes:
+    """Test that different transition types are correctly propagated."""
+
+    def test_different_transition_types(self):
+        """Each branch should use its own transition type for PhaseSpace."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.200, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+
+        # Check that each branch has the correct transition type
+        assert spectrum.branches[0].transition_type == "F2"
+        assert spectrum.branches[1].transition_type == "A"
+
+        # PhaseSpace components should have different W0
+        W, _ = spectrum.get_energy_grid(config)
+        components = spectrum.calculate_components(W)
+
+        branch0_ps = components["branch_0.PhaseSpace"]
+        branch1_ps = components["branch_1.PhaseSpace"]
+
+        # Branch 0 has higher endpoint, so its phase space extends further
+        assert branch0_ps.shape == branch1_ps.shape
+
+    def test_branch_info_transition_type_from_ensdf(self):
+        """BranchInfo should have transition_type field populated from ENSDF."""
+        from beta_spectrum.nuclear_data import BranchInfo
+
+        branch = BranchInfo(
+            level_index=1,
+            level_energy_keV=100.0,
+            intensity=0.5,
+            log_ft=5.0,
+            transition_type="F1",
+        )
+        assert branch.transition_type == "F1"
+        assert branch.level_index == 1
+        assert branch.level_energy_keV == 100.0
+        assert branch.intensity == 0.5
+
+    def test_all_branches_same_transition(self):
+        """All branches can have the same transition type."""
+        config = SpectrumConfig(
+            Z_parent=27,
+            Z_daughter=28,
+            A_number=60,
+            endpoint_MeV=0.320,
+            branches=[
+                BranchConfig(endpoint_MeV=0.310, transition_type="A", intensity=0.5),
+                BranchConfig(endpoint_MeV=0.200, transition_type="A", intensity=0.3),
+                BranchConfig(endpoint_MeV=0.100, transition_type="A", intensity=0.2),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+
+        W, E = spectrum.get_energy_grid(config)
+        values = spectrum(W)
+
+        assert len(values) == len(E)
+        assert np.all(np.isfinite(values))
+        assert np.all(values >= 0)
+
+
+# ---------------------------------------------------------------------------
+# Multi-branch with detector response
+# ---------------------------------------------------------------------------
+
+
+class TestMultiBranchDetectorResponse:
+    """Test multi-branch with detector response convolution."""
+
+    def test_multi_branch_with_detector(self):
+        """Multi-branch spectrum should be convolvable with detector response."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            use_detector_response=True,
+            detector_model="gaussian",
+            detector_sigma_a_keV=1.0,
+            detector_n_channels=256,
+            detector_fano_factor=0.0,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.7),
+                BranchConfig(endpoint_MeV=0.150, transition_type="A", intensity=0.3),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, _ = spectrum.get_energy_grid(config)
+
+        detector = BetaSpectrum.create_detector_from_config(config)
+        convolved = spectrum.convolve_with_detector(detector, W=W, config=config)
+
+        assert convolved.shape == (256,)
+        assert np.all(convolved >= 0), "Convolved spectrum should be non-negative"
+        assert np.sum(convolved) > 0, "Convolved spectrum should have non-zero area"
+
+
+# ---------------------------------------------------------------------------
+# Edge cases
+# ---------------------------------------------------------------------------
+
+
+class TestMultiBranchEdgeCases:
+    """Test edge cases for multi-branch mode."""
+
+    def test_very_low_intensity_branch(self):
+        """Branch with very low intensity should still be computed."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.350,
+            branches=[
+                BranchConfig(endpoint_MeV=0.294, transition_type="F2", intensity=0.999),
+                BranchConfig(endpoint_MeV=0.100, transition_type="A", intensity=0.001),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, _ = spectrum.get_energy_grid(config)
+
+        values = spectrum(W)
+        assert np.all(np.isfinite(values))
+
+    def test_very_high_z_multi_branch(self):
+        """Multi-branch with high Z should be numerically stable."""
+        config = SpectrumConfig(
+            Z_parent=90,
+            Z_daughter=91,
+            A_number=232,
+            endpoint_MeV=1.0,
+            branches=[
+                BranchConfig(endpoint_MeV=0.8, transition_type="A", intensity=0.6),
+                BranchConfig(endpoint_MeV=0.5, transition_type="F1", intensity=0.4),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, _ = spectrum.get_energy_grid(config)
+
+        values = spectrum(W)
+        assert np.all(
+            np.isfinite(values)
+        ), "High-Z multi-branch must not produce NaN/inf"
+
+    def test_no_branches_defaults_to_single(self):
+        """No branches list should default to single-branch mode."""
+        config = SpectrumConfig(
+            Z_parent=43,
+            Z_daughter=44,
+            A_number=99,
+            endpoint_MeV=0.294,
+        )
+        spectrum = BetaSpectrum.from_config(config)
+
+        W, E = spectrum.get_energy_grid(config)
+        values = spectrum(W)
+
+        assert len(values) == len(E)
+        assert np.all(np.isfinite(values))
+
+    def test_branch_values_below_total(self):
+        """Each branch contribution should be <= total spectrum.
+
+        Since total = sum(intensity_i * branch_i) * universal,
+        and all factors are positive, each branch's contribution
+        should not exceed the total.
+        """
+        config = SpectrumConfig(
+            Z_parent=59,
+            Z_daughter=60,
+            A_number=144,
+            endpoint_MeV=2.996,
+            e_step_MeV=0.001,
+            branches=[
+                BranchConfig(
+                    endpoint_MeV=2.996, transition_type="F1", intensity=0.1469
+                ),
+                BranchConfig(
+                    endpoint_MeV=0.341, transition_type="F1", intensity=0.1220
+                ),
+                BranchConfig(
+                    endpoint_MeV=0.253, transition_type="F1", intensity=0.1217
+                ),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, energies = spectrum.get_energy_grid(config)
+
+        total = spectrum(W)
+        bs = spectrum.get_branch_spectra(W)
+        components = spectrum.calculate_components(W)
+
+        # Compute each branch's contribution to the total
+        for i, branch in enumerate(config.branches):
+            universal_product = np.ones_like(W)
+            for name in [
+                "Fermi",
+                "Screening",
+                "FiniteSizeL0",
+                "ChargeDistributionU",
+                "Exchange",
+            ]:
+                if name in components:
+                    universal_product *= components[name]
+            contrib = universal_product * branch.intensity * bs[i]
+
+            # Branch contribution should not exceed total
+            assert np.all(
+                contrib <= total + 1e-10
+            ), f"Branch {i} contribution exceeds total at some points"
+
+            # Branch contribution should be non-negative
+            assert np.all(contrib >= 0), f"Branch {i} contribution has negative values"
+
+    def test_branch_integral_proportional_to_intensity(self):
+        """Integral of each branch contribution should be proportional to its intensity."""
+        config = SpectrumConfig(
+            Z_parent=59,
+            Z_daughter=60,
+            A_number=144,
+            endpoint_MeV=2.996,
+            e_step_MeV=0.001,
+            branches=[
+                BranchConfig(endpoint_MeV=2.996, transition_type="F1", intensity=0.5),
+                BranchConfig(endpoint_MeV=0.5, transition_type="F1", intensity=0.3),
+                BranchConfig(endpoint_MeV=0.3, transition_type="F1", intensity=0.2),
+            ],
+        )
+        spectrum = BetaSpectrum.from_config(config)
+        W, energies = spectrum.get_energy_grid(config)
+
+        bs = spectrum.get_branch_spectra(W)
+        components = spectrum.calculate_components(W)
+
+        integrals = []
+        for i, branch in enumerate(config.branches):
+            universal_product = np.ones_like(W)
+            for name in [
+                "Fermi",
+                "Screening",
+                "FiniteSizeL0",
+                "ChargeDistributionU",
+                "Exchange",
+            ]:
+                if name in components:
+                    universal_product *= components[name]
+            contrib = universal_product * branch.intensity * bs[i]
+            integral = np.trapezoid(contrib, energies)
+            integrals.append(integral)
+
+        total_integral = np.trapezoid(spectrum(W), energies)
+
+        # Each branch integral should be proportional to its intensity
+        # (within numerical tolerance, since different endpoints give different
+        # phase space integrals)
+        for i, (integral, intensity) in enumerate(zip(integrals, config.branches)):
+            assert integral > 0, f"Branch {i} has zero integral"
+
+        # Sum of branch integrals should equal total
+        assert np.isclose(
+            sum(integrals), total_integral, rtol=1e-6
+        ), "Sum of branch integrals should equal total integral"