diff --git a/beta_spectrum/spectrum.py b/beta_spectrum/spectrum.py index 6c107ca..f966b06 100644 --- a/beta_spectrum/spectrum.py +++ b/beta_spectrum/spectrum.py @@ -73,9 +73,9 @@ class SpectrumConfig: e_step_MeV: float = 0.001 neutrino_mass_MeV: float = 0.0 - # Toggle components and corrections - use_phase_space: bool = True - use_fermi: bool = True + # Higher-order corrections. The core spectrum components (phase space + # and the Fermi function, both 0-order) are unconditional: there are + # deliberately no toggles for them. use_screening: bool = True use_finite_size: bool = True use_charge_dist: bool = True @@ -252,26 +252,29 @@ class BetaSpectrum: def _from_config_single( cls, config: SpectrumConfig, logger: Optional[logging.Logger] = None ) -> BetaSpectrum: - """Create single-branch spectrum calculator.""" + """Create single-branch spectrum calculator. + + The core spectrum components (phase space x Fermi function, both + 0-order) are always built; only higher-order corrections are + conditional. + """ components: List[SpectrumComponent] = [] W0 = float(T_to_W(config.endpoint_MeV)) - if config.use_phase_space: - # Convert neutrino mass from MeV to natural units - m_nu_natural = config.neutrino_mass_MeV / ME_MEV - components.append( - PhaseSpace( - W0=W0, - transition_type=config.transition_type, - m_nu=m_nu_natural, - logger=logger, - ) - ) - - if config.use_fermi: - components.append( - FermiFunction(Z=config.Z_daughter, A=config.A_number, logger=logger) + # Core spectrum components (0-order): always included. + # Convert neutrino mass from MeV to natural units + m_nu_natural = config.neutrino_mass_MeV / ME_MEV + components.append( + PhaseSpace( + W0=W0, + transition_type=config.transition_type, + m_nu=m_nu_natural, + logger=logger, ) + ) + components.append( + FermiFunction(Z=config.Z_daughter, A=config.A_number, logger=logger) + ) if config.use_finite_size: components.append( @@ -322,13 +325,13 @@ class BetaSpectrum: Note: Radiative correction is branch-dependent (depends on W0), so it goes into each branch's calculator, not the universal components. """ + if config.branches is None: + raise ValueError("branches must be provided for multi-branch mode") if logger: logger.info("Multi-branch mode: %d branches", len(config.branches)) # Create per-branch calculators branch_spectra: List[BetaSpectrum] = [] - if config.branches is None: - raise ValueError("branches must be provided for multi-branch mode") for i, branch in enumerate(config.branches): if logger: logger.info( @@ -339,37 +342,42 @@ class BetaSpectrum: 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, - neutrino_mass_MeV=config.neutrino_mass_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, - ) + # Per-branch calculator: the per-branch core phase space + # (unconditional) plus optional per-branch radiative + # corrections. The universal Fermi function and universal + # higher-order corrections are applied outside, so they must + # NOT be built here (a Fermi function per branch would + # double-apply it). + W0_branch = float(T_to_W(branch.endpoint_MeV)) + m_nu_natural = config.neutrino_mass_MeV / ME_MEV + branch_components: List[SpectrumComponent] = [ + PhaseSpace( + W0=W0_branch, + transition_type=branch.transition_type, + m_nu=m_nu_natural, + logger=logger, + ) + ] + if config.use_radiative: + branch_components.append( + RadiativeCorrection( + W0=W0_branch, + Z=config.Z_daughter, + A=config.A_number, + use_endpoint_resummation=True, + logger=logger, + ) + ) - branch_spectrum = cls._from_config_single(branch_config, logger) - branch_spectra.append(branch_spectrum) + branch_spectra.append(cls(branch_components, logger=logger)) # 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) - ) + # The Fermi function is a core spectrum component (0-order Coulomb + # correction): always included, not a higher-order correction. + universal_components: List[SpectrumComponent] = [ + FermiFunction(Z=config.Z_daughter, A=config.A_number, logger=logger) + ] if config.use_finite_size: universal_components.append(