From 3216c50943c11cb8ef0d7baba24c33f90d40865e Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Fri, 8 May 2026 08:33:12 +0300 Subject: [PATCH 1/4] feat: per-branch transition types from ENSDF and raw-value plotting - Add transition_type field to BranchInfo and extract per-branch forbiddenness codes from paceENSDF (not defaulting all to 'A') - Fix _multi_branch_call to use np.maximum instead of np.where mask so branches smoothly approach zero at endpoint - Rewrite _plot_multi_branch_spectra and _plot_multi_branch_debug Panel 1 to plot raw (un-normalized) branch contributions so total = sum of branches - Remove independent unit-area normalization that made small branches appear larger than total --- beta_spectrum/nuclear_data.py | 27 +- beta_spectrum/spectrum.py | 695 ++++++++++++++++++++++++++++++++-- 2 files changed, 688 insertions(+), 34 deletions(-) diff --git a/beta_spectrum/nuclear_data.py b/beta_spectrum/nuclear_data.py index 5fb0a0e..c558a9d 100644 --- a/beta_spectrum/nuclear_data.py +++ b/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 - if state[0] == 0 and len(state) > 15 and state[15] is not None: + # Extract per-branch forbiddenness code + 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 - transition_type = FORBIDDENNESS_MAP[code] - # Fallback: use first branch's forbiddenness if no ground state - elif transition_type == "A" and len(state) > 15 and state[15] is not None: - code = str(state[15]) - if code in FORBIDDENNESS_MAP: - forbidden_code = code - transition_type = FORBIDDENNESS_MAP[code] + branch_transition_type = FORBIDDENNESS_MAP[code] + + # Use ground state branch (level_index == 0) for overall forbiddenness + if state[0] == 0: + forbidden_code = ( + str(state[15]) + if len(state) > 15 and state[15] is not None + 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) diff --git a/beta_spectrum/spectrum.py b/beta_spectrum/spectrum.py index fdb4bc9..565e364 100644 --- a/beta_spectrum/spectrum.py +++ b/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] - List of enabled spectral components. + components : List[SpectrumComponent], optional + 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._component_names = [self._get_component_name(c) for c in components] + self.components = components or [] + 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 names - """ + """Get human-readable component name.""" 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) - return components_dict + result[name] = comp(W) + 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( - config.e_step_MeV, config.endpoint_MeV, config.e_step_MeV - ) + if self._is_multi and self.branches: + 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: - self._plot_debug_view(total, components, save_path, commit, timestamp) + if self._is_multi: + 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.""" From 62937100893c949cb392bb5251c350987d7a041a Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Fri, 8 May 2026 08:33:39 +0300 Subject: [PATCH 2/4] feat: multi-branch CLI with intensity cutoff and per-branch info - Add --intensity-cutoff CLI option to filter negligible branches - Fetch DecayInfo from paceENSDF and construct SpectrumConfig with branches including per-branch transition types - Normalize branch intensities to sum to 1.0, filter by cutoff - Display branch info in dry-run output and logging - Add intensity_cutoff to SpectrumConfig --- beta_spectrum/cli.py | 121 ++++++++++++++++++++++++++++++++++++++----- 1 file changed, 109 insertions(+), 12 deletions(-) diff --git a/beta_spectrum/cli.py b/beta_spectrum/cli.py index d566318..19b1e40 100644 --- a/beta_spectrum/cli.py +++ b/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( - "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), - ) + decay_info = get_decay_info_from_paceENSDF(args.nuclide, "beta_minus") + + # Normalize branch intensities to sum to 1.0 + if decay_info.branches: + total_intensity = sum(b.intensity for b in decay_info.branches) + if total_intensity > 0: + branches = [ + BranchConfig( + 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, " From 69eca9773640db06c0aecc3bd7131a4bee166cfa Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Fri, 8 May 2026 08:34:30 +0300 Subject: [PATCH 3/4] test: comprehensive multi-branch decay test suite (30 tests) - BranchConfig validation (intensity bounds, positive endpoint) - Single-branch via multi-branch path compatibility - Multi-branch spectrum calculation and intensity weighting - Per-branch spectra with correct endpoint and transition type - Component extraction (universal vs branch-dependent) - CSV export with total, per-branch, and universal components - Transition type propagation from BranchInfo - Branch values <= total spectrum verification - Branch integral proportionality to intensity - High-Z multi-branch stability - No-branches defaults to single-branch mode --- tests/test_multi_branch.py | 775 +++++++++++++++++++++++++++++++++++++ 1 file changed, 775 insertions(+) create mode 100644 tests/test_multi_branch.py diff --git a/tests/test_multi_branch.py b/tests/test_multi_branch.py new file mode 100644 index 0000000..d9de5a8 --- /dev/null +++ b/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" From ef337073d6e4a248a9186a0ee411b3465f7315d6 Mon Sep 17 00:00:00 2001 From: Evgeniy Unzhakov Date: Fri, 8 May 2026 08:35:01 +0300 Subject: [PATCH 4/4] docs: update TODO and README for B3 completion and BranchConfig export - Cross out B3 in TODO.md with completion details - Export BranchConfig in __init__.py public API - Update README to mention per-branch transition types and intensity cutoff --- README.md | 3 ++- TODO.md | 25 ++++++++++++++++++------- beta_spectrum/__init__.py | 3 ++- 3 files changed, 22 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index e796cf3..c88432e 100644 --- a/README.md +++ b/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 diff --git a/TODO.md b/TODO.md index 2716513..6954b78 100644 --- a/TODO.md +++ b/TODO.md @@ -101,10 +101,21 @@ ______________________________________________________________________ ### B3. Multiple Branch Decays -- [ ] Add support for multiple branch decays - - Calculate total spectral shape weighted by branch intensities - - Handle decays with significant branching to excited states - - Sum contributions from all branches with proper weighting +- [x] ~~Add support for multiple branch decays~~ + - ~~Calculate total spectral shape weighted by branch intensities~~ + - ~~Handle decays with significant branching to excited states~~ + - ~~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\ -**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. +**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. 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). diff --git a/beta_spectrum/__init__.py b/beta_spectrum/__init__.py index 938044a..b836a81 100644 --- a/beta_spectrum/__init__.py +++ b/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",