feat: make phase space and Fermi function unconditional core components

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,11 +252,16 @@ 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:
+        # 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(
@@ -267,8 +272,6 @@ class BetaSpectrum:
                 logger=logger,
             )
         )
-
-        if config.use_fermi:
         components.append(
             FermiFunction(Z=config.Z_daughter, A=config.A_number, logger=logger)
         )
@@ -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,
+            # 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,
-                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,
+                    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(
+        # 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(