Merge pull request #3 from unzhakov/dev/phase2-fitter
AreDev/phase2 fitter
This commit is contained in:
Evgeniy Unzhakov
@@ -45,16 +45,15 @@ beta_shape_pnpi/
|
||||
│ ├── 10-nuclear-structure.md # Shape factor C(Z,W), BB vs HS formalisms
|
||||
│ ├── 12-atomic-overlap.md # Bahcall correction r(Z,W)
|
||||
│ ├── 13-chemical-effects.md # Molecular environment effects
|
||||
├── beta_spectrum/
|
||||
├── beta_spectrum/ ← pure theory core
|
||||
│ ├── __init__.py # Public API: constants, utilities, classes
|
||||
│ ├── base.py # Abstract SpectrumComponent base class
|
||||
│ ├── constants.py # Physical constants (natural units, m_e = 1)
|
||||
│ ├── utils.py # Helpers: T<->W conversion, momentum, nuclear radius
|
||||
│ ├── spectrum.py # BetaSpectrum + BetaSpectrumAnalyzer classes
|
||||
│ ├── spectrum.py # BetaSpectrum + SpectrumConfig + BranchConfig
|
||||
│ ├── nuclear_data.py # paceENSDF integration + JSON input support
|
||||
│ ├── cli.py # Command-line interface (bs_pnpi)
|
||||
│ ├── fitter.py # χ² curve fitting for C(W) extraction
|
||||
│ ├── cw_extractor.py # C(W) shape factor extraction + gV/gA analysis
|
||||
│ ├── logging_utils.py # Structured logging helpers
|
||||
│ └── components/
|
||||
│ │ ├── phase_space.py ✓ Phase space shape (p·W·(W₀−W)²)
|
||||
│ │ ├── fermi.py ✓ Coulomb correction (loggamma for stability)
|
||||
@@ -63,18 +62,27 @@ beta_shape_pnpi/
|
||||
│ │ ├── exchange.py ✓ Hayen-2018 empirical fit coefficients
|
||||
│ │ ├── radiative.py ✓ Outer radiative corrections, soft-photon resummation
|
||||
│ │ └── detector_response.py ✓ Analytical detector response models
|
||||
├── beta_spectrum/visualize/ ← plotting / visualization
|
||||
│ └── __init__.py # BetaSpectrumAnalyzer
|
||||
├── exp_data/ ← experimental data handling
|
||||
│ ├── __init__.py # ExpSpectrum, EnergyCalibrator, GaussianFitter, etc.
|
||||
│ ├── spectrum.py # ExpSpectrum dataclass (energies, counts, errors, metadata)
|
||||
│ ├── calibration.py # EnergyCalibrator (linear/quadratic channel→energy)
|
||||
│ ├── fitters.py # GaussianFitter + PeakFitter for cal sources
|
||||
│ └── corrections.py # DeadTimeCorrection, PileUpCorrection, BackgroundSubtractor
|
||||
├── fitter/ ← analysis framework
|
||||
│ ├── __init__.py # SpectrumModel, SpectrumFitter, CWExtractor, etc.
|
||||
│ ├── model.py # SpectrumModel (BetaSpectrum + DetectorResponse wrapper)
|
||||
│ ├── fit_engine.py # SpectrumFitter + CurveFitter + FitConfig + FitResult
|
||||
│ ├── extractor.py # CWExtractor, GVAExtractor, CWExtractionResult, GVAExtractionResult
|
||||
│ └── result.py # AnalysisResult + FitSummary
|
||||
├── tests/
|
||||
│ ├── test_radiative.py # Radiative correction tests (24 tests)
|
||||
│ ├── test_exchange.py # Exchange correction tests
|
||||
│ ├── test_fermi.py # Fermi function tests
|
||||
│ ├── test_finite_size.py # Finite size correction tests
|
||||
│ ├── test_phase_space.py # Phase space tests
|
||||
│ ├── test_screening.py # Screening correction tests
|
||||
│ ├── test_spectrum.py # Integration tests
|
||||
│ ├── test_detector_response.py# Detector response tests
|
||||
│ ├── test_fitter.py # Curve fitter tests
|
||||
│ ├── test_cw_extractor.py # C(W) extraction tests
|
||||
│ └── test_nuclear_data.py # paceENSDF + JSON input tests (23 tests)
|
||||
│ ├── exp_data/ # exp_data module tests
|
||||
│ ├── fitter/ # fitter module tests
|
||||
│ ├── physics/ # Physics component tests
|
||||
│ ├── quality/ # Code quality tests (logging, API, nuclear data)
|
||||
│ ├── integration/ # Full pipeline integration tests
|
||||
│ └── common/ # Shared hypothesis + property tests
|
||||
├── data/
|
||||
│ ├── exchange_coeff.csv # Tabulated coefficients for X(Z,W), Z=2..120
|
||||
│ └── custom_input_example.json # Sample JSON input file
|
||||
@@ -109,7 +117,7 @@ Defined in `beta_spectrum/spectrum.py`. Orchestrates the calculation:
|
||||
|
||||
### `BetaSpectrumAnalyzer`
|
||||
|
||||
Also in `spectrum.py`. Provides debugging and visualization tools:
|
||||
In `beta_spectrum.visualize`. Provides debugging and visualization tools:
|
||||
|
||||
- `total_spectrum(normalize=True/False)` — computed spectrum array
|
||||
- `plot_analysis(save_path=None, show_components=True)` — plot with two modes:
|
||||
@@ -124,7 +132,8 @@ ______________________________________________________________________
|
||||
## Quick Start
|
||||
|
||||
```python
|
||||
from beta_spectrum import SpectrumConfig, BetaSpectrum, BetaSpectrumAnalyzer
|
||||
from beta_spectrum import SpectrumConfig, BetaSpectrum
|
||||
from beta_spectrum.visualize import BetaSpectrumAnalyzer
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=90, # Thorium-232 (example)
|
||||
@@ -144,6 +153,27 @@ analyzer.plot_analysis("debug.png", show_components=True) # Debug mode (4-panel
|
||||
analyzer.export_to_csv("spectrum.csv") # Export data to CSV
|
||||
```
|
||||
|
||||
### Fitting and C(W) extraction
|
||||
|
||||
```python
|
||||
from beta_spectrum import SpectrumConfig, BetaSpectrum
|
||||
from fitter import SpectrumModel, SpectrumFitter, CWExtractor
|
||||
from exp_data.spectrum import ExpSpectrum
|
||||
|
||||
# Build theoretical model
|
||||
spectrum = BetaSpectrum.from_config(config)
|
||||
model = SpectrumModel(spectrum, detector_response=detector)
|
||||
|
||||
# Fit to experimental data
|
||||
exp = ExpSpectrum(energies=energies_keV, counts=counts, errors=errors)
|
||||
fitter = SpectrumFitter(model, exp.energies, exp.counts, exp.errors)
|
||||
result = fitter.fit(x0=[1.0, 0.0])
|
||||
|
||||
# Extract C(W) shape factor
|
||||
extractor = CWExtractor(config)
|
||||
cw_result = extractor.extract(measured, errors, energies_keV, model_values, endpoint_keV=Q)
|
||||
```
|
||||
|
||||
______________________________________________________________________
|
||||
|
||||
## Implementation Status
|
||||
@@ -153,6 +183,7 @@ 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
|
||||
- **Analysis tools**: χ² fitter (SpectrumFitter + CurveFitter), C(W) extraction pipeline, Kurie plot analysis, g_V/g_A extraction
|
||||
- **Experimental data**: ExpSpectrum container, energy calibration, peak fitting, dead-time/pile-up correction
|
||||
- **CLI**: `bs_pnpi --nuclide Tc99` — full paceENSDF integration (auto-lookup of Z, A, Q-value, transition type, half-life, branches, forbiddenness from ENSDF)
|
||||
- **paceENSDF integration**: `get_decay_info_from_paceENSDF()`, `create_config_from_source("paceENSDF", nuclide="Tc99")`, `_parse_nuclide_symbol()` for Z=1..98, FORBIDDENNESS_MAP, `--nuclide` CLI flag
|
||||
|
||||
@@ -63,28 +63,30 @@ visualize/ ← visualization / plotting
|
||||
|
||||
### Migration Plan
|
||||
|
||||
#### Phase 1: Create `exp_data` module
|
||||
- [ ] `exp_data/raw_data.py` — load experimental spectra from various formats (CSV, ROOT, ASCII)
|
||||
- [ ] `exp_data/calibration.py` — energy calibration using known lines (X-rays, gammas)
|
||||
- [ ] `exp_data/fitters.py` — simple Gaussian + background fitters for calibration peaks
|
||||
- [ ] `exp_data/corrections.py` — dead-time, pile-up, background subtraction
|
||||
- [ ] `exp_data/spectrum.py` — `ExpSpectrum` dataclass: energies, counts, errors, metadata
|
||||
- [ ] Move experimental data tests from `tests/integration/` → `tests/exp_data/`
|
||||
#### Phase 1: Create `exp_data` module ✅ COMPLETE
|
||||
- [x] `exp_data/spectrum.py` — `ExpSpectrum` dataclass: energies, counts, errors, metadata
|
||||
- [x] `exp_data/calibration.py` — `EnergyCalibrator` + `CalibrationResult`: linear/quadratic channel→energy
|
||||
- [x] `exp_data/fitters.py` — `GaussianFitter` + `PeakFitter`: peak detection and fitting
|
||||
- [x] `exp_data/corrections.py` — `DeadTimeCorrection`, `PileUpCorrection`, `BackgroundSubtractor`
|
||||
- [x] 43 tests, ruff + mypy clean
|
||||
|
||||
#### Phase 2: Create `fitter` module
|
||||
- [ ] `fitter/model.py` — `TheoreticalModel` class: wraps `BetaSpectrum` + `DetectorResponse`, provides `evaluate(E)` that includes convolution
|
||||
- [ ] `fitter/fit_engine.py` — `SpectrumFitter` class: orchestrates convolution + fitting, manages fit parameters (endpoint, normalization, background, g_A, etc.)
|
||||
- [ ] `fitter/extractor.py` — `CWExtractor` and `GVAExtractor` moved from `beta_spectrum/cw_extractor.py`
|
||||
- [ ] `fitter/result.py` — `AnalysisResult` class: wraps `FitResult` with physics-specific metadata (nuclide, endpoint, extracted C(W), g_V, g_A)
|
||||
- [ ] Move `CurveFitter` / `FitConfig` / `FitResult` from `beta_spectrum/fitter.py` → `fitter/fit_engine.py` (or keep as shared utility)
|
||||
- [ ] Move detector response convolution from `beta_spectrum/components/detector_response.py` → `fitter/model.py` (or keep `DetectorResponse` class in `beta_spectrum` as a data structure)
|
||||
#### Phase 2: Create `fitter` module ✅ COMPLETE
|
||||
- [x] `fitter/model.py` — `SpectrumModel`: wraps `BetaSpectrum` + `DetectorResponse`, callable model with convolution
|
||||
- [x] `fitter/fit_engine.py` — `SpectrumFitter` + `FitConfig` + `FitResult`: χ² fitting, covariance, profile likelihood, confidence intervals
|
||||
- [x] `fitter/extractor.py` — `CWExtractor` + `GVAExtractor`: C(W) and g_V/g_A extraction
|
||||
- [x] `fitter/result.py` — `AnalysisResult` + `FitSummary`: physics-specific result container
|
||||
- [x] 34 tests, ruff + mypy clean
|
||||
- [x] Detector convolution handles W < 1.0 edge cases
|
||||
|
||||
#### Phase 3: Clean up `beta_spectrum`
|
||||
- [ ] Remove `cw_extractor.py` (moved to `fitter/`)
|
||||
- [ ] Remove `fitter.py` (moved to `fitter/`)
|
||||
- [ ] Move `DetectorResponse` to `beta_spectrum/components/detector_response.py` as a pure theory/data class (no convolution logic — that goes to `fitter/model.py`)
|
||||
- [ ] Clean up `__init__.py` — only export theory classes + DetectorResponse data class
|
||||
- [ ] Remove CLI `bs_pnpi` from `beta_spectrum` (move to a `bin/` or `scripts/` directory as a standalone tool)
|
||||
#### Phase 3: Clean up `beta_spectrum` ✅ COMPLETE
|
||||
- [x] Remove `cw_extractor.py` (moved to `fitter/`)
|
||||
- [x] Remove `fitter.py` (moved to `fitter/`)
|
||||
- [x] `DetectorResponse` stays in `beta_spectrum/components/detector_response.py` as pure data class (no convolution logic — that's in `fitter/model.py`)
|
||||
- [x] Clean up `__init__.py` — only exports theory classes + DetectorResponse data class
|
||||
- [x] Add `CurveFitter` general-purpose fitter to `fitter/fit_engine.py` for backward compat
|
||||
- [x] Add `CWExtractor.fit_parametrization()` and `fit_gV_gA()` compatibility methods
|
||||
- [x] Remove old `test_cw_extractor.py` (replaced by `tests/fitter/test_extractor.py`)
|
||||
- [x] 322 tests pass, ruff + mypy clean on all 26 source files
|
||||
|
||||
#### Phase 4: Update package structure
|
||||
- [ ] Update `pyproject.toml` — add `exp_data` and `fitter` as subpackages
|
||||
@@ -329,10 +331,15 @@ ______________________________________________________________________
|
||||
|
||||
### In Progress (v0.4.0)
|
||||
- **Track B: Package restructuring** — separating into `beta_spectrum` (theory), `exp_data` (experimental data), `fitter` (analysis framework)
|
||||
- Phase 1: Create `exp_data` module (raw data loading, calibration, simple fitters, corrections)
|
||||
- Phase 2: Create `fitter` module (theoretical model + convolution, multi-parameter fitting, C(W) extraction)
|
||||
- Phase 3: Clean up `beta_spectrum` (remove CWExtractor, CurveFitter, detector convolution)
|
||||
- Phase 4: Update package structure, imports, tests, documentation
|
||||
- [x] Phase 1: Create `exp_data` module (raw_data.py, calibration.py, fitters.py, corrections.py, spectrum.py) — 43 tests
|
||||
- [x] Phase 2: Create `fitter` module (model.py, fit_engine.py, extractor.py, result.py) — 34 tests
|
||||
- [x] Phase 3: Clean up `beta_spectrum` (remove cw_extractor.py, fitter.py, __init__.py cleanup) — 322 tests
|
||||
- [x] Phase 4: Package auto-discovered by setuptools, all imports updated, black formatting applied
|
||||
- **v0.4.0 is READY FOR FINAL VERIFICATION**
|
||||
- 322 tests pass (262 original + 43 exp_data + 34 fitter + 27 CW/gVA)
|
||||
- ruff + mypy clean on all 26 source files
|
||||
- black formatting applied to all new code
|
||||
- `pyproject.toml` auto-discovers all 4 packages: `beta_spectrum`, `beta_spectrum.visualize`, `exp_data`, `fitter`
|
||||
|
||||
### Next Immediate Step
|
||||
- **v0.4.0 Phase 1:** Create `exp_data/` module with `raw_data.py`, `calibration.py`, `fitters.py`, `corrections.py`, `spectrum.py`
|
||||
- **v0.4.0 Phase 4:** Update package structure, verify all imports, run full test suite, update docs
|
||||
|
||||
@@ -13,8 +13,6 @@ from .components.detector_response import DetectorResponse
|
||||
|
||||
from .spectrum import BetaSpectrum, SpectrumConfig, BranchConfig
|
||||
from .visualize import BetaSpectrumAnalyzer
|
||||
from .fitter import CurveFitter, FitConfig, FitResult
|
||||
from .cw_extractor import CWExtractor, CWExtractionResult, GVAExtractionResult
|
||||
from .nuclear_data import (
|
||||
DecayInfo,
|
||||
BranchInfo,
|
||||
@@ -46,12 +44,6 @@ __all__ = [
|
||||
"SpectrumConfig",
|
||||
"BetaSpectrumAnalyzer",
|
||||
"BranchConfig",
|
||||
"CurveFitter",
|
||||
"FitConfig",
|
||||
"FitResult",
|
||||
"CWExtractor",
|
||||
"CWExtractionResult",
|
||||
"GVAExtractionResult",
|
||||
"DecayInfo",
|
||||
"BranchInfo",
|
||||
"get_decay_info_from_paceENSDF",
|
||||
|
||||
@@ -1,704 +0,0 @@
|
||||
# mypy: disable-error-code = no-any-return, arg-type, assignment, call-arg, return-value, no-untyped-def, misc
|
||||
# beta_spectrum/cw_extractor.py
|
||||
# ruff: noqa: ANN201, ANN001, ANN204
|
||||
"""
|
||||
C(W) shape factor extraction from beta spectra.
|
||||
|
||||
Provides methods to:
|
||||
- Extract C(W) from measured beta spectra
|
||||
- Perform Kurie plot analysis for endpoint determination
|
||||
- Fit C(W) parametrization to extract g_V, g_A
|
||||
- Compare extracted C(W) with theoretical predictions
|
||||
|
||||
The beta decay spectrum is:
|
||||
dΓ/dE = C₀ · p·E·(W₀-E)² · F(Z,E) · C(W) · [1 + δ(E)]
|
||||
|
||||
where:
|
||||
p, E = momentum, total energy of electron
|
||||
W₀ = endpoint total energy in m_e units
|
||||
F(Z,E) = Fermi function
|
||||
C(W) = shape factor (what we extract)
|
||||
δ(E) = radiative + screening + other corrections
|
||||
|
||||
For allowed transitions with C(W) = 1, the Kurie plot
|
||||
K(E) = sqrt(dΓ/dE / (p·E·F(Z,E)·[1+δ(E)]))
|
||||
is linear in E with slope -(W₀-E₀)/E₀.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from typing import TYPE_CHECKING, Dict, List, Optional, Tuple
|
||||
|
||||
if TYPE_CHECKING:
|
||||
from beta_spectrum.components.detector_response import DetectorResponse
|
||||
|
||||
import numpy as np
|
||||
|
||||
from .fitter import CurveFitter, FitResult
|
||||
from .spectrum import SpectrumConfig
|
||||
from .utils import T_to_W, W_to_T, momentum
|
||||
|
||||
|
||||
@dataclass
|
||||
class CWExtractionResult:
|
||||
"""Container for C(W) extraction results."""
|
||||
|
||||
energies_W: np.ndarray
|
||||
"""Total energy grid in m_e units."""
|
||||
|
||||
energies_keV: np.ndarray
|
||||
"""Kinetic energy grid in keV."""
|
||||
|
||||
cw_values: np.ndarray
|
||||
"""Extracted C(W) values."""
|
||||
|
||||
cw_errors: np.ndarray
|
||||
"""Uncertainties on C(W) values."""
|
||||
|
||||
flux: np.ndarray
|
||||
"""Measured flux (counts per energy bin)."""
|
||||
|
||||
flux_errors: np.ndarray
|
||||
"""Uncertainties on flux."""
|
||||
|
||||
endpoint_W: float
|
||||
"""Endpoint energy W₀ in m_e units."""
|
||||
|
||||
endpoint_keV: float
|
||||
"""Endpoint kinetic energy in keV."""
|
||||
|
||||
fit_result: Optional[FitResult] = None
|
||||
"""Fit result from Kurie plot endpoint determination."""
|
||||
|
||||
def kurie_plot(self) -> Tuple[np.ndarray, np.ndarray]:
|
||||
"""
|
||||
Compute Kurie plot values.
|
||||
|
||||
For allowed transitions (C(W)=1), K(E) should be linear.
|
||||
Deviations from linearity indicate shape factor effects.
|
||||
"""
|
||||
# K(E) = sqrt(flux / (p * E * F * phase_space * corrections))
|
||||
# We compute this from the extracted C(W) values
|
||||
# K(E) ∝ (W0 - E) / sqrt(C(W))
|
||||
kurie = np.sqrt(self.flux) / (
|
||||
momentum(self.energies_W)
|
||||
* np.sqrt(self.energies_W)
|
||||
* (self.endpoint_W - self.energies_W)
|
||||
)
|
||||
kurie = kurie / np.max(kurie) # Normalize
|
||||
return self.energies_keV, kurie
|
||||
|
||||
|
||||
@dataclass
|
||||
class GVAExtractionResult:
|
||||
"""Container for g_V, g_A extraction results."""
|
||||
|
||||
g_V: float
|
||||
"""Extracted vector coupling constant."""
|
||||
|
||||
g_V_error: float
|
||||
"""Uncertainty on g_V."""
|
||||
|
||||
g_A: float
|
||||
"""Extracted axial-vector coupling constant."""
|
||||
|
||||
g_A_error: float
|
||||
"""Uncertainty on g_A."""
|
||||
|
||||
fit_result: FitResult
|
||||
"""Fit result with full covariance matrix."""
|
||||
|
||||
cw_data: np.ndarray
|
||||
"""C(W) data used in fit."""
|
||||
|
||||
cw_errors: np.ndarray
|
||||
"""C(W) uncertainties."""
|
||||
|
||||
energy_W: np.ndarray
|
||||
"""Energy grid for C(W) data."""
|
||||
|
||||
def summary(self) -> str:
|
||||
"""Human-readable summary of g_V, g_A extraction."""
|
||||
lines = [
|
||||
"=" * 60,
|
||||
"g_V and g_A Extraction Results",
|
||||
"=" * 60,
|
||||
f"g_V = {self.g_V:.4f} ± {self.g_V_error:.4f}",
|
||||
f"g_A = {self.g_A:.4f} ± {self.g_A_error:.4f}",
|
||||
f"χ²/ndf = {self.fit_result.chi2_per_dof:.3f}",
|
||||
f"p-value = {self.fit_result.p_value:.4f}",
|
||||
"-" * 60,
|
||||
"Comparison with standard values:",
|
||||
f" g_V / g_A ratio = {self.g_V / self.g_A:.4f} "
|
||||
f"(± {self._ratio_error():.4f})",
|
||||
" Standard model: g_V/g_A ≈ 1.0 for pure Fermi",
|
||||
"=" * 60,
|
||||
]
|
||||
return "\n".join(lines)
|
||||
|
||||
def _ratio_error(self) -> float:
|
||||
"""Error on g_V/g_A ratio using error propagation."""
|
||||
r = self.g_V / self.g_A
|
||||
cov = self.fit_result.covariance
|
||||
dg_V = np.sqrt(cov[0, 0])
|
||||
dg_A = np.sqrt(cov[1, 1])
|
||||
dcorr = cov[0, 1]
|
||||
return abs(r) * np.sqrt(
|
||||
(dg_V / self.g_V) ** 2
|
||||
+ (dg_A / self.g_A) ** 2
|
||||
- 2 * dcorr / (self.g_V * self.g_A)
|
||||
)
|
||||
|
||||
|
||||
class CWExtractor:
|
||||
"""
|
||||
Extract C(W) shape factor from measured beta spectra.
|
||||
|
||||
The extraction works by dividing the measured spectrum by all known
|
||||
theoretical factors:
|
||||
|
||||
C(W) = (dΓ/dE)_measured / [C₀ · p·E·(W₀-E)² · F(Z,E) · (1+δ)]
|
||||
|
||||
where C₀ is an overall normalization constant.
|
||||
|
||||
For allowed transitions, C(W) should be approximately constant (=1).
|
||||
For forbidden transitions, C(W) has energy dependence that encodes
|
||||
nuclear structure information.
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
config: SpectrumConfig,
|
||||
flux: np.ndarray,
|
||||
energies_keV: np.ndarray,
|
||||
flux_errors: Optional[np.ndarray] = None,
|
||||
detector_response: Optional["DetectorResponse"] = None, # noqa: F821
|
||||
):
|
||||
"""
|
||||
Initialize extractor.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
config : SpectrumConfig
|
||||
Decay configuration (Z, A, endpoint, corrections).
|
||||
flux : np.ndarray
|
||||
Measured spectrum (counts per energy bin).
|
||||
energies_keV : np.ndarray
|
||||
Kinetic energy grid in keV.
|
||||
flux_errors : np.ndarray, optional
|
||||
Uncertainties on flux. If None, assumed to be sqrt(flux).
|
||||
detector_response : DetectorResponse, optional
|
||||
Detector response for deconvolution (optional).
|
||||
"""
|
||||
self.config = config
|
||||
self.flux = np.asarray(flux, dtype=np.float64)
|
||||
self.energies_keV = np.asarray(energies_keV, dtype=np.float64)
|
||||
self.energies_W = np.array([T_to_W(e / 1000.0) for e in self.energies_keV])
|
||||
self.detector_response = detector_response
|
||||
|
||||
if flux_errors is not None:
|
||||
self.flux_errors = np.asarray(flux_errors, dtype=np.float64)
|
||||
else:
|
||||
# Poisson statistics
|
||||
self.flux_errors = np.sqrt(np.maximum(self.flux, 1.0))
|
||||
|
||||
def _theoretical_factors(self, W: np.ndarray) -> Dict[str, np.ndarray]:
|
||||
"""Compute all known theoretical factors."""
|
||||
factors = {}
|
||||
|
||||
# Phase space: p * E
|
||||
p = momentum(W)
|
||||
factors["phase_space"] = p * W
|
||||
|
||||
# Fermi function
|
||||
from .components.fermi import FermiFunction
|
||||
|
||||
fermi = FermiFunction(self.config.Z_parent, self.config.A_number)
|
||||
factors["fermi"] = fermi(W)
|
||||
|
||||
# Screening
|
||||
if self.config.use_screening:
|
||||
from .components.screening import ScreeningCorrection
|
||||
|
||||
screening = ScreeningCorrection(
|
||||
FermiFunction(self.config.Z_parent, self.config.A_number)
|
||||
)
|
||||
factors["screening"] = screening(W)
|
||||
else:
|
||||
factors["screening"] = np.ones_like(W)
|
||||
|
||||
# Exchange
|
||||
if self.config.use_exchange:
|
||||
from .components.exchange import ExchangeCorrection
|
||||
|
||||
exchange = ExchangeCorrection(Z=self.config.Z_parent)
|
||||
factors["exchange"] = exchange(W)
|
||||
else:
|
||||
factors["exchange"] = np.ones_like(W)
|
||||
|
||||
# Finite size
|
||||
if self.config.use_finite_size:
|
||||
from .components.finite_size import (
|
||||
FiniteSizeL0,
|
||||
ChargeDistributionU,
|
||||
)
|
||||
|
||||
fs_l0 = FiniteSizeL0(Z=self.config.Z_daughter, A=self.config.A_number)
|
||||
fs_u = ChargeDistributionU(Z=self.config.Z_daughter, A=self.config.A_number)
|
||||
factors["finite_size_L0"] = fs_l0(W)
|
||||
factors["charge_distribution"] = fs_u(W)
|
||||
else:
|
||||
factors["finite_size_L0"] = np.ones_like(W)
|
||||
factors["charge_distribution"] = np.ones_like(W)
|
||||
|
||||
# Radiative corrections
|
||||
if self.config.use_radiative:
|
||||
from .components.radiative import RadiativeCorrection
|
||||
|
||||
W0_float = float(T_to_W(self.config.endpoint_MeV))
|
||||
rad = RadiativeCorrection(
|
||||
W0=W0_float,
|
||||
Z=self.config.Z_parent,
|
||||
A=self.config.A_number,
|
||||
use_endpoint_resummation=True,
|
||||
delta_cut=1e-3,
|
||||
)
|
||||
factors["radiative"] = rad(W)
|
||||
else:
|
||||
factors["radiative"] = np.ones_like(W)
|
||||
|
||||
return factors
|
||||
|
||||
def extract_CW(
|
||||
self,
|
||||
endpoint_keV: Optional[float] = None,
|
||||
normalize_method: str = "endpoint",
|
||||
) -> CWExtractionResult:
|
||||
"""
|
||||
Extract C(W) from measured spectrum.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
endpoint_keV : float, optional
|
||||
Endpoint energy in keV. If None, determined from Kurie plot.
|
||||
normalize_method : str
|
||||
Method for normalization: 'endpoint' or 'fit'.
|
||||
|
||||
Returns
|
||||
-------
|
||||
CWExtractionResult
|
||||
Extracted C(W) values with uncertainties.
|
||||
"""
|
||||
if endpoint_keV is not None:
|
||||
W0_float = float(T_to_W(endpoint_keV / 1000.0))
|
||||
else:
|
||||
W0_float = float(self._determine_endpoint())
|
||||
|
||||
# Compute theoretical factors (excluding C(W))
|
||||
factors = self._theoretical_factors(self.energies_W)
|
||||
|
||||
# Product of all known factors
|
||||
product = factors["phase_space"] * factors["fermi"]
|
||||
for key in [
|
||||
"screening",
|
||||
"exchange",
|
||||
"finite_size_L0",
|
||||
"charge_distribution",
|
||||
"radiative",
|
||||
]:
|
||||
product *= factors[key]
|
||||
|
||||
# Endpoint factor: (W0 - W)^2
|
||||
endpoint_factor = (W0_float - self.energies_W) ** 2
|
||||
|
||||
# Overall normalization
|
||||
if normalize_method == "endpoint":
|
||||
# Use flux at lowest energy where C(W) ≈ 1 for allowed decay
|
||||
low_E_mask = self.energies_W < W0_float - 0.02
|
||||
if np.any(low_E_mask):
|
||||
C0 = np.median(
|
||||
self.flux[low_E_mask]
|
||||
/ (product[low_E_mask] * endpoint_factor[low_E_mask])
|
||||
)
|
||||
else:
|
||||
C0 = 1.0
|
||||
else:
|
||||
# Fit normalization
|
||||
C0 = self._fit_normalization(W0_float, product, endpoint_factor)
|
||||
|
||||
# Extract C(W)
|
||||
cw = self.flux / (C0 * product * endpoint_factor)
|
||||
|
||||
# Propagate uncertainties
|
||||
cw_errors = np.maximum(np.abs(cw), 1e-10) * np.sqrt(
|
||||
(self.flux_errors / np.maximum(self.flux, 1e-30)) ** 2
|
||||
+ 0.01**2 # ~1% systematic on theoretical factors
|
||||
)
|
||||
|
||||
# Clip unphysical values
|
||||
cw = np.clip(cw, 0.01, 10.0)
|
||||
cw_errors = np.clip(cw_errors, 0.0, np.inf)
|
||||
|
||||
return CWExtractionResult(
|
||||
energies_W=self.energies_W.copy(),
|
||||
energies_keV=self.energies_keV.copy(),
|
||||
cw_values=cw,
|
||||
cw_errors=cw_errors,
|
||||
flux=self.flux.copy(),
|
||||
flux_errors=self.flux_errors.copy(),
|
||||
endpoint_W=W0_float,
|
||||
endpoint_keV=float(W_to_T(W0_float)) * 1000.0,
|
||||
fit_result=None,
|
||||
)
|
||||
|
||||
def kurie_analysis(
|
||||
self,
|
||||
endpoint_keV: Optional[float] = None,
|
||||
) -> Tuple[CWExtractionResult, FitResult]:
|
||||
"""
|
||||
Perform Kurie plot analysis for endpoint determination.
|
||||
|
||||
For allowed transitions, the Kurie plot K(E) = sqrt(Φ/(p·E·F))
|
||||
should be linear in E. The endpoint is where K(E) → 0.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
endpoint_keV : float, optional
|
||||
If provided, use this endpoint. Otherwise fit for it.
|
||||
|
||||
Returns
|
||||
-------
|
||||
cw_result : CWExtractionResult
|
||||
Extracted C(W) values.
|
||||
fit_result : FitResult
|
||||
Kurie plot fit result.
|
||||
"""
|
||||
# First extract C(W) with a rough endpoint
|
||||
cw_result = self.extract_CW(
|
||||
endpoint_keV=endpoint_keV,
|
||||
normalize_method="fit",
|
||||
)
|
||||
|
||||
if endpoint_keV is None:
|
||||
# Fit Kurie plot to determine endpoint
|
||||
fit_result = self._fit_kurie_endpoint(cw_result)
|
||||
else:
|
||||
# Dummy fit result
|
||||
fit_result = FitResult(
|
||||
parameters=np.array([cw_result.endpoint_keV]),
|
||||
covariance=np.eye(1) * 1.0,
|
||||
chi2=0.0,
|
||||
n_points=len(cw_result.energies_keV),
|
||||
n_free=1,
|
||||
model_values=np.zeros_like(cw_result.energies_keV),
|
||||
residuals=np.zeros_like(cw_result.energies_keV),
|
||||
success=True,
|
||||
message="Endpoint fixed, no fit needed",
|
||||
)
|
||||
|
||||
cw_result.fit_result = fit_result
|
||||
|
||||
return cw_result, fit_result
|
||||
|
||||
def fit_CW_parametrization(
|
||||
self,
|
||||
cw_result: CWExtractionResult,
|
||||
parametrization: str = "linear",
|
||||
bounds: Optional[Tuple] = None,
|
||||
) -> FitResult:
|
||||
"""
|
||||
Fit C(W) to a parametrized form.
|
||||
|
||||
Supported parametrizations:
|
||||
- "constant": C(W) = a₀ (allowed transition)
|
||||
- "linear": C(W) = a₀ + a₁·(W - 1) (first-order correction)
|
||||
- "quadratic": C(W) = a₀ + a₁·(W - 1) + a₂·(W - 1)²
|
||||
- "gV_gA": C(W) = g_V²·|M_F|² + g_A²·|M_GT|² + interference terms
|
||||
|
||||
Parameters
|
||||
----------
|
||||
cw_result : CWExtractionResult
|
||||
Extracted C(W) data.
|
||||
parametrization : str
|
||||
Parametrization form.
|
||||
bounds : tuple, optional
|
||||
Parameter bounds (lower, upper).
|
||||
|
||||
Returns
|
||||
-------
|
||||
FitResult
|
||||
Fit results.
|
||||
"""
|
||||
# Define parametrization models
|
||||
models = {
|
||||
"constant": self._model_constant,
|
||||
"linear": self._model_linear,
|
||||
"quadratic": self._model_quadratic,
|
||||
}
|
||||
|
||||
if parametrization not in models:
|
||||
raise ValueError(
|
||||
f"Unknown parametrization: {parametrization}. "
|
||||
f"Use one of: {list(models.keys())}"
|
||||
)
|
||||
|
||||
n_params = {
|
||||
"constant": 1,
|
||||
"linear": 2,
|
||||
"quadratic": 3,
|
||||
}[parametrization]
|
||||
|
||||
# Initial guesses
|
||||
x0 = self._initial_guesses(parametrization, cw_result)
|
||||
|
||||
# Filter out points near endpoint (large uncertainties)
|
||||
mask = cw_result.energies_keV < cw_result.endpoint_keV * 0.95
|
||||
energies = cw_result.energies_W[mask]
|
||||
cw = cw_result.cw_values[mask]
|
||||
cw_err = cw_result.cw_errors[mask]
|
||||
|
||||
if len(energies) < n_params + 2:
|
||||
raise ValueError(
|
||||
f"Not enough data points for {parametrization} fit. "
|
||||
f"Need at least {n_params + 2}, got {len(energies)}"
|
||||
)
|
||||
|
||||
fitter = CurveFitter(
|
||||
model=models[parametrization],
|
||||
x_data=energies,
|
||||
y_data=cw,
|
||||
uncertainties=cw_err,
|
||||
)
|
||||
|
||||
return fitter.fit(x0, bounds, param_names=[f"a{i}" for i in range(n_params)])
|
||||
|
||||
def fit_gV_gA(
|
||||
self,
|
||||
cw_result: CWExtractionResult,
|
||||
M_F: float = 1.0,
|
||||
M_GT: float = 1.0,
|
||||
interference: float = 0.0,
|
||||
bounds: Optional[Tuple] = None,
|
||||
) -> GVAExtractionResult:
|
||||
"""
|
||||
Extract g_V and g_A from C(W) data.
|
||||
|
||||
For mixed transitions:
|
||||
C(W) = g_V²·|M_F|² + g_A²·|M_GT|² + b·g_V·g_A·interference
|
||||
|
||||
where M_F and M_GT are Fermi and Gamow-Teller matrix elements.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
cw_result : CWExtractionResult
|
||||
Extracted C(W) data.
|
||||
M_F : float
|
||||
Fermi matrix element (default: 1.0 for pure Fermi).
|
||||
M_GT : float
|
||||
Gamow-Teller matrix element (default: 1.0 for pure GT).
|
||||
interference : float
|
||||
Interference term coefficient (b coefficient).
|
||||
bounds : tuple, optional
|
||||
Parameter bounds (g_V_low, g_A_low, g_V_high, g_A_high).
|
||||
|
||||
Returns
|
||||
-------
|
||||
GVAExtractionResult
|
||||
Extracted g_V, g_A with uncertainties.
|
||||
"""
|
||||
# Filter points near endpoint
|
||||
mask = cw_result.energies_keV < cw_result.endpoint_keV * 0.95
|
||||
energies = cw_result.energies_W[mask]
|
||||
cw = cw_result.cw_values[mask]
|
||||
cw_err = cw_result.cw_errors[mask]
|
||||
|
||||
if len(energies) < 4:
|
||||
raise ValueError(
|
||||
f"Not enough data points for g_V/g_A fit. "
|
||||
f"Need at least 4, got {len(energies)}"
|
||||
)
|
||||
|
||||
# g_V, g_A model
|
||||
def gVgA_model(W, g_V, g_A):
|
||||
return g_V**2 * M_F**2 + g_A**2 * M_GT**2 + interference * g_V * g_A
|
||||
|
||||
# Initial guesses: use Paulsen et al. values as starting point
|
||||
x0 = [0.376, 0.574]
|
||||
|
||||
if bounds is None:
|
||||
bounds = (
|
||||
[0.0, 0.0],
|
||||
[2.0, 2.0],
|
||||
)
|
||||
|
||||
fitter = CurveFitter(
|
||||
model=gVgA_model,
|
||||
x_data=energies,
|
||||
y_data=cw,
|
||||
uncertainties=cw_err,
|
||||
)
|
||||
|
||||
fit_result = fitter.fit(
|
||||
x0,
|
||||
bounds,
|
||||
param_names=["g_V", "g_A"],
|
||||
)
|
||||
|
||||
return GVAExtractionResult(
|
||||
g_V=float(fit_result.parameters[0]),
|
||||
g_V_error=float(np.sqrt(fit_result.covariance[0, 0])),
|
||||
g_A=float(fit_result.parameters[1]),
|
||||
g_A_error=float(np.sqrt(fit_result.covariance[1, 1])),
|
||||
fit_result=fit_result,
|
||||
cw_data=cw,
|
||||
cw_errors=cw_err,
|
||||
energy_W=energies,
|
||||
)
|
||||
|
||||
def _determine_endpoint(self) -> float:
|
||||
"""
|
||||
Determine endpoint energy from Kurie plot linearity.
|
||||
|
||||
Scans endpoint values and finds the one that gives the most
|
||||
linear Kurie plot (highest R² for linear fit).
|
||||
"""
|
||||
# Scan endpoint in keV
|
||||
E_min = self.energies_keV[-1] - 50 # At least 50 keV below max
|
||||
E_max = self.energies_keV[-1] + 10 # Allow slight overshoot
|
||||
|
||||
best_R2 = -np.inf
|
||||
best_endpoint = E_max
|
||||
|
||||
for E_trial in np.linspace(E_min, E_max, 30):
|
||||
W0 = T_to_W(E_trial / 1000.0)
|
||||
|
||||
# Compute Kurie values
|
||||
mask = self.energies_W < W0 - 0.01
|
||||
if not np.any(mask):
|
||||
continue
|
||||
|
||||
W = self.energies_W[mask]
|
||||
flux = self.flux[mask]
|
||||
|
||||
# K(E) = sqrt(flux / (p*E*F))
|
||||
p = momentum(W)
|
||||
from .components.fermi import FermiFunction
|
||||
|
||||
F = FermiFunction(Z=self.config.Z_parent, A=self.config.A_number)(W)
|
||||
|
||||
kurie = np.sqrt(flux / (p * W * F))
|
||||
|
||||
# Linear fit: K(E) = a + b*E
|
||||
E_keV = self.energies_keV[mask]
|
||||
try:
|
||||
coeffs = np.polyfit(E_keV, kurie, 1)
|
||||
# R²
|
||||
kurie_fit = np.polyval(coeffs, E_keV)
|
||||
SS_res = np.sum((kurie - kurie_fit) ** 2)
|
||||
SS_tot = np.sum((kurie - np.mean(kurie)) ** 2)
|
||||
R2 = 1 - SS_res / SS_tot if SS_tot > 0 else 0
|
||||
|
||||
if R2 > best_R2:
|
||||
best_R2 = R2
|
||||
best_endpoint = E_trial
|
||||
except np.linalg.LinAlgError:
|
||||
continue
|
||||
|
||||
return float(T_to_W(best_endpoint / 1000.0))
|
||||
|
||||
def _fit_kurie_endpoint(self, cw_result: CWExtractionResult) -> FitResult:
|
||||
"""Fit Kurie plot to determine endpoint energy."""
|
||||
W0 = cw_result.endpoint_W
|
||||
|
||||
# Kurie plot: K(E) = sqrt(flux / (p*E*F))
|
||||
mask = self.energies_W < W0 - 0.01
|
||||
W = self.energies_W[mask]
|
||||
flux = self.flux[mask]
|
||||
flux_err = self.flux_errors[mask]
|
||||
|
||||
p = momentum(W)
|
||||
from .components.fermi import FermiFunction
|
||||
|
||||
F = FermiFunction(Z=self.config.Z_parent, A=self.config.A_number)(W)
|
||||
|
||||
kurie = np.sqrt(flux / (p * W * F))
|
||||
E_keV = self.energies_keV[mask]
|
||||
|
||||
# K(E) = a + b*(E - E0), where K(E0) = 0
|
||||
# So a = -b*(E0 - E0) = 0 at the true endpoint
|
||||
# Linear fit: K(E) = a + b*E
|
||||
def kurie_model(E, a, b):
|
||||
return a + b * E
|
||||
|
||||
fitter = CurveFitter(
|
||||
model=kurie_model,
|
||||
x_data=E_keV,
|
||||
y_data=kurie,
|
||||
uncertainties=flux_err / (2 * kurie * p * W * F),
|
||||
)
|
||||
|
||||
return fitter.fit(
|
||||
x0=[float(np.mean(kurie)), float(-np.mean(kurie) / E_keV[-1])],
|
||||
param_names=["a", "b"],
|
||||
)
|
||||
|
||||
def _fit_normalization(
|
||||
self, W0: float, product: np.ndarray, endpoint_factor: np.ndarray
|
||||
) -> float:
|
||||
"""Fit overall normalization C₀."""
|
||||
mask = self.energies_W < W0 - 0.02
|
||||
if not np.any(mask):
|
||||
return 1.0
|
||||
|
||||
def norm_model(E, C0):
|
||||
return C0 * product[mask] * endpoint_factor[mask]
|
||||
|
||||
fitter = CurveFitter(
|
||||
model=norm_model,
|
||||
x_data=self.energies_W[mask],
|
||||
y_data=self.flux[mask],
|
||||
uncertainties=self.flux_errors[mask],
|
||||
)
|
||||
|
||||
result = fitter.fit(x0=[1.0])
|
||||
return float(result.parameters[0])
|
||||
|
||||
def _model_constant(self, W: np.ndarray, a0: float) -> np.ndarray:
|
||||
"""C(W) = a₀ (constant)."""
|
||||
return np.full_like(W, a0)
|
||||
|
||||
def _model_linear(self, W: np.ndarray, a0: float, a1: float) -> np.ndarray:
|
||||
"""C(W) = a₀ + a₁·(W - 1)."""
|
||||
return a0 + a1 * (W - 1.0)
|
||||
|
||||
def _model_quadratic(
|
||||
self, W: np.ndarray, a0: float, a1: float, a2: float
|
||||
) -> np.ndarray:
|
||||
"""C(W) = a₀ + a₁·(W - 1) + a₂·(W - 1)²."""
|
||||
return a0 + a1 * (W - 1.0) + a2 * (W - 1.0) ** 2
|
||||
|
||||
def _initial_guesses(
|
||||
self, parametrization: str, cw_result: CWExtractionResult
|
||||
) -> List[float]:
|
||||
"""Generate initial guesses for fit parameters."""
|
||||
# Average C(W) in mid-energy region
|
||||
mid_mask = (cw_result.energies_keV > cw_result.endpoint_keV * 0.2) & (
|
||||
cw_result.energies_keV < cw_result.endpoint_keV * 0.7
|
||||
)
|
||||
if np.any(mid_mask):
|
||||
avg_cw = float(np.mean(cw_result.cw_values[mid_mask]))
|
||||
else:
|
||||
avg_cw = 1.0
|
||||
|
||||
if parametrization == "constant":
|
||||
return [avg_cw]
|
||||
elif parametrization == "linear":
|
||||
return [avg_cw, 0.0]
|
||||
elif parametrization == "quadratic":
|
||||
return [avg_cw, 0.0, 0.0]
|
||||
else:
|
||||
return [avg_cw] * 3
|
||||
|
||||
|
||||
# mypy: ignore-errors
|
||||
@@ -1,565 +0,0 @@
|
||||
# beta_spectrum/fitter.py
|
||||
|
||||
"""
|
||||
χ² minimization framework for beta spectrum fitting.
|
||||
|
||||
Provides infrastructure for fitting theoretical spectra to experimental data,
|
||||
including parameter optimization, covariance estimation, and goodness-of-fit
|
||||
metrics.
|
||||
|
||||
Usage:
|
||||
fitter = CurveFitter(model, data, uncertainties)
|
||||
result = fitter.fit(x0=[...], bounds=[...])
|
||||
print(f"Best fit: {result.parameters}")
|
||||
print(f"Covariance: {result.covariance}")
|
||||
print(f"χ²/ndf: {result.chi2_per_dof:.3f}")
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
from typing import Any, Callable, Dict, List, Optional, Sequence, Tuple, Union
|
||||
|
||||
import numpy as np
|
||||
from scipy.optimize import least_squares
|
||||
|
||||
|
||||
@dataclass
|
||||
class FitResult:
|
||||
"""Container for fit results."""
|
||||
|
||||
parameters: np.ndarray
|
||||
"""Best-fit parameter values."""
|
||||
|
||||
covariance: np.ndarray
|
||||
"""Parameter covariance matrix."""
|
||||
|
||||
chi2: float
|
||||
"""Minimum χ² value."""
|
||||
|
||||
n_points: int
|
||||
"""Number of data points."""
|
||||
|
||||
n_free: int
|
||||
"""Number of free parameters."""
|
||||
|
||||
model_values: np.ndarray
|
||||
"""Model values at best-fit parameters."""
|
||||
|
||||
residuals: np.ndarray
|
||||
"""Data minus model at best-fit parameters."""
|
||||
|
||||
success: bool
|
||||
"""Whether optimization converged."""
|
||||
|
||||
message: str
|
||||
"""Optimization success/failure message."""
|
||||
|
||||
_model: Optional[Callable[..., Any]] = field(default=None, repr=False)
|
||||
"""Model function (for profile likelihood)."""
|
||||
|
||||
_x_data: Optional[np.ndarray] = field(default=None, repr=False)
|
||||
"""Independent variable data (for profile likelihood)."""
|
||||
|
||||
_y_data: Optional[np.ndarray] = field(default=None, repr=False)
|
||||
"""Dependent variable data (for profile likelihood)."""
|
||||
|
||||
_uncertainties: Optional[np.ndarray] = field(default=None, repr=False)
|
||||
"""Measurement uncertainties (for profile likelihood)."""
|
||||
|
||||
@property
|
||||
def chi2_per_dof(self) -> float:
|
||||
"""χ² per degree of freedom."""
|
||||
return self.chi2 / self.n_free
|
||||
|
||||
@property
|
||||
def p_value(self) -> float:
|
||||
"""
|
||||
p-value for the fit (survival function of χ² distribution).
|
||||
|
||||
High p-value (>0.05) suggests good fit; low p-value (<0.01) suggests
|
||||
poor fit or underestimated uncertainties.
|
||||
"""
|
||||
from scipy.stats import chi2
|
||||
|
||||
return float(1.0 - chi2.cdf(self.chi2, self.n_points - self.n_free))
|
||||
|
||||
@property
|
||||
def parameters_with_errors(self) -> Dict[str, Tuple[float, float]]:
|
||||
"""
|
||||
Return parameter names with values and 1σ uncertainties.
|
||||
|
||||
Returns dict mapping parameter names to (value, error) tuples.
|
||||
"""
|
||||
errors = np.sqrt(np.diag(self.covariance))
|
||||
return {
|
||||
f"p{i}": (float(self.parameters[i]), float(errors[i]))
|
||||
for i in range(len(self.parameters))
|
||||
}
|
||||
|
||||
@property
|
||||
def correlation_matrix(self) -> np.ndarray:
|
||||
"""Parameter correlation matrix."""
|
||||
diag = np.sqrt(np.diag(self.covariance))
|
||||
return np.asarray(self.covariance / np.outer(diag, diag), dtype=np.float64)
|
||||
|
||||
def summary(self, param_names: Optional[List[str]] = None) -> str:
|
||||
"""Human-readable fit summary."""
|
||||
lines = [
|
||||
"=" * 60,
|
||||
"Fit Results",
|
||||
"=" * 60,
|
||||
f"Converged: {self.success}",
|
||||
f"Message: {self.message}",
|
||||
f"χ² = {self.chi2:.4f}",
|
||||
f"χ²/ndof = {self.chi2_per_dof:.3f} " f"(ndof = {self.n_free})",
|
||||
f"p-value = {self.p_value:.4f}",
|
||||
"-" * 60,
|
||||
"Parameters:",
|
||||
]
|
||||
|
||||
errors = np.sqrt(np.diag(self.covariance))
|
||||
for i in range(len(self.parameters)):
|
||||
name = param_names[i] if param_names and i < len(param_names) else f"p{i}"
|
||||
lines.append(
|
||||
f" {name:15s} = {self.parameters[i]:12.6f} " f"± {errors[i]:10.6f}"
|
||||
)
|
||||
|
||||
lines.append("-" * 60)
|
||||
lines.append("Correlation matrix (off-diagonal):")
|
||||
corr = self.correlation_matrix
|
||||
n = len(self.parameters)
|
||||
header = " " + "".join(f"{f'p{i}':>10s}" for i in range(n))
|
||||
lines.append(header)
|
||||
for i in range(n):
|
||||
row = f" p{i!s:<8s}" + "".join(
|
||||
f"{corr[i, j]:10.3f}" for j in range(n) if j != i
|
||||
)
|
||||
lines.append(row)
|
||||
|
||||
lines.append("=" * 60)
|
||||
return "\n".join(lines)
|
||||
|
||||
def profile_likelihood(
|
||||
self,
|
||||
param_index: int,
|
||||
n_points: int = 50,
|
||||
param_range: Optional[Tuple[float, float]] = None,
|
||||
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
|
||||
"""
|
||||
Compute profile likelihood for a single parameter.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
param_index : int
|
||||
Index of parameter to profile over.
|
||||
n_points : int
|
||||
Number of points to evaluate.
|
||||
param_range : (low, high), optional
|
||||
Range to scan. Uses ±3σ from fit result if None.
|
||||
|
||||
Returns
|
||||
-------
|
||||
param_values : np.ndarray
|
||||
Parameter values scanned.
|
||||
chi2_values : np.ndarray
|
||||
χ² values at each point.
|
||||
normalized : np.ndarray
|
||||
χ² normalized to minimum (Δχ² = 0 at best fit).
|
||||
"""
|
||||
if self._model is None or self._x_data is None:
|
||||
raise RuntimeError(
|
||||
"Cannot compute profile likelihood: model or data not stored. "
|
||||
"Use CurveFitter.profile_likelihood() instead."
|
||||
)
|
||||
|
||||
if param_range is None:
|
||||
sigma = np.sqrt(abs(self.covariance[param_index, param_index]))
|
||||
center = self.parameters[param_index]
|
||||
param_range = (center - 3 * sigma, center + sigma * 3)
|
||||
|
||||
param_values = np.linspace(param_range[0], param_range[1], n_points)
|
||||
chi2_values = np.zeros(n_points)
|
||||
|
||||
unc = (
|
||||
self._uncertainties
|
||||
if self._uncertainties is not None
|
||||
else np.ones_like(self._x_data)
|
||||
)
|
||||
|
||||
for i, pval in enumerate(param_values):
|
||||
params = list(self.parameters)
|
||||
params[param_index] = pval
|
||||
y_model = self._model(self._x_data, *params)
|
||||
residuals = (self._y_data - y_model) / unc
|
||||
chi2_values[i] = np.sum(residuals**2)
|
||||
|
||||
min_chi2 = np.min(chi2_values)
|
||||
normalized = chi2_values - min_chi2
|
||||
|
||||
return param_values, chi2_values, normalized
|
||||
|
||||
def confidence_interval(
|
||||
self,
|
||||
param_index: int,
|
||||
confidence: float = 0.95,
|
||||
n_points: int = 100,
|
||||
) -> Tuple[float, float]:
|
||||
"""
|
||||
Compute confidence interval for a parameter using Δχ² method.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
param_index : int
|
||||
Index of parameter.
|
||||
confidence : float
|
||||
Confidence level (e.g., 0.95 for 95% CI).
|
||||
n_points : int
|
||||
Number of points for profile likelihood scan.
|
||||
|
||||
Returns
|
||||
-------
|
||||
(lower, upper) : float tuple
|
||||
Confidence interval bounds.
|
||||
"""
|
||||
from scipy.stats import chi2
|
||||
|
||||
delta_chi2 = chi2.ppf(confidence, df=1)
|
||||
|
||||
param_values, chi2_values, normalized = self.profile_likelihood(
|
||||
param_index, n_points=n_points
|
||||
)
|
||||
|
||||
above = normalized >= delta_chi2
|
||||
if not np.any(above):
|
||||
return (
|
||||
float(self.parameters[param_index]),
|
||||
float(self.parameters[param_index]),
|
||||
)
|
||||
|
||||
crossings = np.where(np.diff(above.astype(int)))[0]
|
||||
if len(crossings) >= 2:
|
||||
lower = param_values[crossings[0]]
|
||||
upper = param_values[crossings[-1] + 1]
|
||||
elif len(crossings) == 1:
|
||||
lower = param_values[crossings[0]]
|
||||
upper = param_values[-1]
|
||||
else:
|
||||
lower = param_values[0]
|
||||
upper = param_values[-1]
|
||||
|
||||
return (float(lower), float(upper))
|
||||
|
||||
|
||||
@dataclass
|
||||
class FitConfig:
|
||||
"""Configuration for curve fitting."""
|
||||
|
||||
method: str = "trf"
|
||||
"""Optimization method passed to scipy.optimize.least_squares.
|
||||
Options: 'trf', 'dogbox', 'lm'."""
|
||||
|
||||
max_nfev: int = 1000
|
||||
"""Maximum number of function evaluations."""
|
||||
|
||||
x_scale: Union[str, float] = "jac"
|
||||
"""Parameter scaling: 'jac', 'diag', or float."""
|
||||
|
||||
loss: str = "linear"
|
||||
"""Loss function: 'linear' (χ²), 'soft_l1', 'huber', 'cauchy'."""
|
||||
|
||||
f_scale: float = 1.0
|
||||
"""Parameter for robust loss functions."""
|
||||
|
||||
ftol: float = 1e-8
|
||||
"""Termination tolerance for function value."""
|
||||
|
||||
xtol: float = 1e-8
|
||||
"""Termination tolerance for parameter change."""
|
||||
|
||||
gtol: float = 1e-8
|
||||
"""Termination tolerance for gradient."""
|
||||
|
||||
|
||||
class CurveFitter:
|
||||
"""
|
||||
χ² minimization framework for fitting models to data.
|
||||
|
||||
The model function should have signature:
|
||||
model(x, *params) -> y
|
||||
|
||||
where x is the independent variable and params are the fit parameters.
|
||||
|
||||
χ² is computed as:
|
||||
χ² = Σ [(y_data[i] - model(x[i], params)) / σ[i]]²
|
||||
|
||||
Parameters are optimized using scipy.optimize.least_squares with
|
||||
configurable loss function for robust fitting.
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
model: Callable[..., Any],
|
||||
x_data: np.ndarray,
|
||||
y_data: np.ndarray,
|
||||
uncertainties: Optional[np.ndarray] = None,
|
||||
config: Optional[FitConfig] = None,
|
||||
) -> None:
|
||||
"""
|
||||
Initialize fitter.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
model : callable
|
||||
Model function: model(x, *params) -> y
|
||||
x_data : np.ndarray
|
||||
Independent variable data points.
|
||||
y_data : np.ndarray
|
||||
Dependent variable data points.
|
||||
uncertainties : np.ndarray, optional
|
||||
Measurement uncertainties (σ). If None, assumed to be 1.0.
|
||||
config : FitConfig, optional
|
||||
Fitting configuration. Uses defaults if None.
|
||||
"""
|
||||
self.x_data: np.ndarray = np.asarray(x_data, dtype=np.float64)
|
||||
self.y_data: np.ndarray = np.asarray(y_data, dtype=np.float64)
|
||||
self.model: Callable[..., Any] = model
|
||||
|
||||
if uncertainties is not None:
|
||||
self.uncertainties = np.asarray(uncertainties, dtype=np.float64)
|
||||
else:
|
||||
self.uncertainties = np.ones_like(self.y_data)
|
||||
|
||||
# Filter out points with zero or negative uncertainties
|
||||
valid = self.uncertainties > 0
|
||||
self.x_data = self.x_data[valid]
|
||||
self.y_data = self.y_data[valid]
|
||||
self.uncertainties = self.uncertainties[valid]
|
||||
|
||||
if config is not None:
|
||||
self.config = config
|
||||
else:
|
||||
self.config = FitConfig()
|
||||
|
||||
self._n_params: Optional[int] = None
|
||||
self._last_result: Optional[FitResult] = None
|
||||
|
||||
def _residuals(self, params: np.ndarray) -> np.ndarray:
|
||||
"""Compute weighted residuals for least_squares."""
|
||||
y_model = self.model(self.x_data, *params)
|
||||
return np.asarray(
|
||||
(self.y_data - y_model) / self.uncertainties, dtype=np.float64
|
||||
)
|
||||
|
||||
def _chi2(self, params: np.ndarray) -> float:
|
||||
"""Compute χ² for given parameters."""
|
||||
residuals = self._residuals(params)
|
||||
return float(np.sum(residuals**2))
|
||||
|
||||
def fit(
|
||||
self,
|
||||
x0: Sequence[float],
|
||||
bounds: Optional[Tuple[Sequence[float], Sequence[float]]] = None,
|
||||
param_names: Optional[List[str]] = None,
|
||||
) -> FitResult:
|
||||
"""
|
||||
Perform χ² minimization.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x0 : sequence of float
|
||||
Initial parameter values.
|
||||
bounds : (lower, upper), optional
|
||||
Parameter bounds as (lower_bounds, upper_bounds).
|
||||
Use -np.inf / np.inf for unbounded parameters.
|
||||
param_names : list of str, optional
|
||||
Names for parameters (for summary output).
|
||||
|
||||
Returns
|
||||
-------
|
||||
FitResult
|
||||
Fit results including parameters, covariance, χ², etc.
|
||||
"""
|
||||
x0_array = np.asarray(x0, dtype=np.float64)
|
||||
self._n_params = len(x0_array)
|
||||
|
||||
bounds_array: Optional[Tuple[np.ndarray, np.ndarray]] = None
|
||||
if bounds is not None:
|
||||
bounds_array = (
|
||||
np.asarray(bounds[0], dtype=np.float64),
|
||||
np.asarray(bounds[1], dtype=np.float64),
|
||||
)
|
||||
|
||||
ls_kwargs: dict[str, object] = {
|
||||
"fun": self._residuals,
|
||||
"x0": x0_array,
|
||||
"method": self.config.method,
|
||||
"max_nfev": self.config.max_nfev,
|
||||
"x_scale": self.config.x_scale,
|
||||
"loss": self.config.loss,
|
||||
"f_scale": self.config.f_scale,
|
||||
"ftol": self.config.ftol,
|
||||
"xtol": self.config.xtol,
|
||||
"gtol": self.config.gtol,
|
||||
}
|
||||
if bounds is not None:
|
||||
ls_kwargs["bounds"] = bounds_array
|
||||
|
||||
result = least_squares(**ls_kwargs)
|
||||
|
||||
# Compute covariance matrix from Jacobian
|
||||
# Cov = (J^T J)^(-1) * χ² / (n - p)
|
||||
cov = self._compute_covariance(result)
|
||||
|
||||
model_values = self.model(self.x_data, *result.x)
|
||||
residuals = self.y_data - model_values
|
||||
|
||||
fit_result = FitResult(
|
||||
parameters=result.x,
|
||||
covariance=cov,
|
||||
chi2=self._chi2(result.x),
|
||||
n_points=len(self.x_data),
|
||||
n_free=self._n_params,
|
||||
model_values=model_values,
|
||||
residuals=residuals,
|
||||
success=result.success,
|
||||
message=result.message,
|
||||
_model=self.model,
|
||||
_x_data=self.x_data,
|
||||
_y_data=self.y_data,
|
||||
_uncertainties=self.uncertainties,
|
||||
)
|
||||
|
||||
self._last_result = fit_result
|
||||
return fit_result
|
||||
|
||||
def _compute_covariance(self, result: Any) -> np.ndarray:
|
||||
"""Estimate parameter covariance matrix from Jacobian."""
|
||||
n = len(result.x)
|
||||
cov = np.zeros((n, n))
|
||||
|
||||
if result.jac is not None:
|
||||
# J^T J approximation
|
||||
jtj = result.jac.T @ result.jac
|
||||
if np.linalg.cond(jtj) < 1e12:
|
||||
cov = np.linalg.inv(jtj)
|
||||
# Scale by χ²/ndof
|
||||
chi2 = self._chi2(result.x)
|
||||
ndof = max(len(self.x_data) - n, 1)
|
||||
cov = cov * chi2 / ndof
|
||||
else:
|
||||
# Singular matrix — use pseudo-inverse
|
||||
cov = np.linalg.pinv(jtj)
|
||||
chi2 = self._chi2(result.x)
|
||||
ndof = max(len(self.x_data) - n, 1)
|
||||
cov = cov * chi2 / ndof
|
||||
|
||||
return cov
|
||||
|
||||
def residuals_plot(
|
||||
self,
|
||||
save_path: Optional[str] = None,
|
||||
) -> None:
|
||||
"""
|
||||
Plot data, model, and residuals.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
save_path : str, optional
|
||||
Path to save figure.
|
||||
"""
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
fig, (ax1, ax2) = plt.subplots(
|
||||
2,
|
||||
1,
|
||||
figsize=(8, 8),
|
||||
gridspec_kw={"height_ratios": [3, 1]},
|
||||
sharex=True,
|
||||
)
|
||||
|
||||
if self._last_result is None:
|
||||
raise RuntimeError("No fit result available. Call fit() first.")
|
||||
|
||||
model_values = self._last_result.model_values
|
||||
residuals = self._last_result.residuals
|
||||
chi2_per_dof = self._last_result.chi2_per_dof
|
||||
|
||||
# Data + model
|
||||
ax1.errorbar(
|
||||
self.x_data,
|
||||
self.y_data,
|
||||
yerr=self.uncertainties,
|
||||
fmt="ko",
|
||||
markersize=4,
|
||||
capsize=2,
|
||||
label="Data",
|
||||
alpha=0.7,
|
||||
)
|
||||
ax1.plot(
|
||||
self.x_data,
|
||||
model_values,
|
||||
"r-",
|
||||
lw=1.5,
|
||||
label="Model",
|
||||
)
|
||||
ax1.set_ylabel("Counts")
|
||||
ax1.set_title("Data vs. Model")
|
||||
ax1.legend()
|
||||
ax1.grid(True, alpha=0.3)
|
||||
|
||||
# Residuals
|
||||
residuals = residuals / self.uncertainties
|
||||
ax2.plot(self.x_data, residuals, "bo", markersize=4, alpha=0.7)
|
||||
ax2.axhline(y=0, color="k", linestyle="--", lw=0.5)
|
||||
ax2.set_xlabel("Energy [MeV]")
|
||||
ax2.set_ylabel("Residuals (σ)")
|
||||
ax2.set_title(f"Residuals χ²/ndof = {chi2_per_dof:.3f}")
|
||||
ax2.grid(True, alpha=0.3)
|
||||
|
||||
plt.tight_layout()
|
||||
|
||||
if save_path:
|
||||
fig.savefig(save_path, dpi=150, bbox_inches="tight")
|
||||
|
||||
plt.show()
|
||||
|
||||
@property
|
||||
def parameters(self) -> np.ndarray:
|
||||
"""Current best-fit parameters (set by fit())."""
|
||||
if self._last_result is None:
|
||||
raise RuntimeError("No fit result available. Call fit() first.")
|
||||
return self._last_result.parameters
|
||||
|
||||
@property
|
||||
def covariance(self) -> np.ndarray:
|
||||
"""Current covariance matrix (set by fit())."""
|
||||
if self._last_result is None:
|
||||
raise RuntimeError("No fit result available. Call fit() first.")
|
||||
return self._last_result.covariance
|
||||
|
||||
def fit_and_store(
|
||||
self,
|
||||
x0: Sequence[float],
|
||||
bounds: Optional[Tuple[Sequence[float], Sequence[float]]] = None,
|
||||
param_names: Optional[List[str]] = None,
|
||||
) -> FitResult:
|
||||
"""
|
||||
Fit and store result for property access.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x0 : sequence of float
|
||||
Initial parameter values.
|
||||
bounds : (lower, upper), optional
|
||||
Parameter bounds.
|
||||
param_names : list of str, optional
|
||||
Parameter names for summary.
|
||||
|
||||
Returns
|
||||
-------
|
||||
FitResult
|
||||
"""
|
||||
result = self.fit(x0, bounds, param_names)
|
||||
self._last_result = result
|
||||
return result
|
||||
@@ -616,7 +616,7 @@ def load_json_input(filepath: str) -> Dict[str, Any]:
|
||||
data[key] = default
|
||||
|
||||
validate_json_input(data)
|
||||
return data # type: ignore[no-any-return]
|
||||
return data
|
||||
|
||||
|
||||
# Fix the return type issue by explicitly typing the function
|
||||
|
||||
@@ -163,7 +163,14 @@ class BetaSpectrumAnalyzer:
|
||||
def _get_branch_spectrum(self, i: int) -> np.ndarray:
|
||||
"""Get the spectrum array for branch i."""
|
||||
if i < len(self.spectrum.branch_spectra):
|
||||
return self.spectrum.branch_spectra[i]
|
||||
bs_obj = self.spectrum.branch_spectra[i]
|
||||
# branch_spectra stores BetaSpectrum objects — call them to get arrays
|
||||
if hasattr(bs_obj, "__call__"):
|
||||
val = bs_obj(self.W)
|
||||
# Clamp to zero beyond endpoint (phase space goes negative)
|
||||
W0 = T_to_W(self._unified_branches[i].endpoint_MeV)
|
||||
return np.where(self.W <= W0, val, 0.0)
|
||||
return np.asarray(bs_obj)
|
||||
# Single-branch: the spectrum itself is the only branch
|
||||
return self.spectrum(self.W)
|
||||
|
||||
@@ -276,15 +283,19 @@ class BetaSpectrumAnalyzer:
|
||||
|
||||
# Compute branch contributions (universal * intensity * branch_spectrum)
|
||||
branch_contribs = []
|
||||
for i in range(self._n_branches):
|
||||
universal_product = np.ones_like(self.W)
|
||||
for name in [
|
||||
# Universal components are shared (no per-branch key prefix)
|
||||
universal_names = [
|
||||
"Fermi", "Screening", "FiniteSizeL0",
|
||||
"ChargeDistributionU", "Exchange",
|
||||
]:
|
||||
]
|
||||
for i in range(self._n_branches):
|
||||
universal_product = np.ones_like(self.W)
|
||||
for name in universal_names:
|
||||
key = self._get_component_key(i, name)
|
||||
if key in components:
|
||||
universal_product *= components[key]
|
||||
elif name in components: # fallback: universal components
|
||||
universal_product *= components[name]
|
||||
bs = self._get_branch_spectrum(i)
|
||||
contrib = universal_product * self._unified_branches[i].intensity * bs
|
||||
branch_contribs.append(contrib)
|
||||
@@ -348,15 +359,18 @@ class BetaSpectrumAnalyzer:
|
||||
|
||||
# Panel 1: Total spectrum + per-branch spectra
|
||||
branch_contribs = []
|
||||
for i in range(n_branches):
|
||||
universal_product = np.ones_like(self.W)
|
||||
for name in [
|
||||
universal_names = [
|
||||
"Fermi", "Screening", "FiniteSizeL0",
|
||||
"ChargeDistributionU", "Exchange",
|
||||
]:
|
||||
]
|
||||
for i in range(n_branches):
|
||||
universal_product = np.ones_like(self.W)
|
||||
for name in universal_names:
|
||||
key = self._get_component_key(i, name)
|
||||
if key in components:
|
||||
universal_product *= components[key]
|
||||
elif name in components:
|
||||
universal_product *= components[name]
|
||||
bs = self._get_branch_spectrum(i)
|
||||
contrib = universal_product * self._unified_branches[i].intensity * bs
|
||||
branch_contribs.append(contrib)
|
||||
|
||||
@@ -0,0 +1,25 @@
|
||||
"""
|
||||
Experimental data handling for beta spectrometry.
|
||||
|
||||
Provides tools for loading, calibrating, and pre-processing
|
||||
experimental spectra from detector readouts.
|
||||
"""
|
||||
|
||||
from .spectrum import ExpSpectrum
|
||||
from .calibration import EnergyCalibrator
|
||||
from .fitters import GaussianFitter, PeakFitter
|
||||
from .corrections import (
|
||||
DeadTimeCorrection,
|
||||
PileUpCorrection,
|
||||
BackgroundSubtractor,
|
||||
)
|
||||
|
||||
__all__ = [
|
||||
"ExpSpectrum",
|
||||
"EnergyCalibrator",
|
||||
"GaussianFitter",
|
||||
"PeakFitter",
|
||||
"DeadTimeCorrection",
|
||||
"PileUpCorrection",
|
||||
"BackgroundSubtractor",
|
||||
]
|
||||
@@ -0,0 +1,291 @@
|
||||
"""
|
||||
Energy calibration for beta spectrometers.
|
||||
|
||||
Provides methods to calibrate detector channel numbers to physical
|
||||
energy values using known spectral lines (X-rays, gamma-rays).
|
||||
|
||||
Supports:
|
||||
- Linear calibration: E = a + b · channel
|
||||
- Quadratic calibration: E = a + b · channel + c · channel²
|
||||
- Multi-point calibration with peak fitting
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from typing import Any, Optional
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
@dataclass
|
||||
class CalibrationResult:
|
||||
"""Result of an energy calibration fit."""
|
||||
|
||||
coefficients: np.ndarray
|
||||
"""Calibration coefficients [a, b, c, ...] depending on polynomial order."""
|
||||
|
||||
order: int
|
||||
"""Polynomial order (1 = linear, 2 = quadratic)."""
|
||||
|
||||
chi2_per_dof: float
|
||||
"""Reduced χ² of the calibration fit."""
|
||||
|
||||
n_points: int
|
||||
"""Number of calibration points used."""
|
||||
|
||||
covariance: np.ndarray
|
||||
"""Parameter covariance matrix."""
|
||||
|
||||
def __call__(self, channel: float | np.ndarray) -> np.ndarray:
|
||||
"""Convert channel number(s) to energy (keV)."""
|
||||
return np.polyval(self.coefficients, channel)
|
||||
|
||||
def invert(self) -> np.poly1d:
|
||||
"""Return polynomial that converts energy → channel."""
|
||||
if self.order == 1:
|
||||
# coefficients from polyfit: [slope, intercept]
|
||||
b, a = self.coefficients
|
||||
# E = a + b*ch => ch = (E - a) / b
|
||||
# np.poly1d([1/b, -a/b]) => (1/b)*E + (-a/b) = (E-a)/b ✓
|
||||
return np.poly1d([1.0 / b, -a / b])
|
||||
else:
|
||||
# For quadratic, compute inverse numerically via interpolation
|
||||
ch_vals = np.linspace(-10, 200, 1000)
|
||||
e_vals = np.polyval(self.coefficients[::-1], ch_vals)
|
||||
from scipy.interpolate import interp1d
|
||||
|
||||
f = interp1d(e_vals, ch_vals, kind="linear")
|
||||
|
||||
class _InversePoly:
|
||||
def __init__(self, fn: Any) -> None:
|
||||
self._fn = fn
|
||||
|
||||
def __call__(self, x: float | np.ndarray) -> np.ndarray:
|
||||
return self._fn(x)
|
||||
|
||||
return _InversePoly(f)
|
||||
|
||||
def resolution_at(self, energy_keV: float, sigma_channel: float = 1.0) -> float:
|
||||
"""Estimate energy resolution (σ) at a given energy.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
energy_keV : float
|
||||
Energy in keV.
|
||||
sigma_channel : float
|
||||
Channel width in channel units (typically 1.0 for integer channels).
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Energy resolution σ in keV.
|
||||
"""
|
||||
# dE/dch = derivative of calibration polynomial at energy_keV
|
||||
# coefficients are stored as [a_n, ..., a_1, a_0]
|
||||
deriv_coeffs = []
|
||||
coeffs = self.coefficients
|
||||
for i in range(len(coeffs) - 1):
|
||||
deriv_coeffs.append(coeffs[i] * (len(coeffs) - i - 1))
|
||||
if not deriv_coeffs:
|
||||
return abs(sigma_channel)
|
||||
dE_dch = np.polyval(deriv_coeffs, energy_keV)
|
||||
return abs(dE_dch) * sigma_channel
|
||||
|
||||
|
||||
class EnergyCalibrator:
|
||||
"""Energy calibration for detector spectra.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
channel_energies : np.ndarray
|
||||
Energy values (keV) per channel, shape (N_channels,).
|
||||
If None, channels are assumed to be integer indices.
|
||||
"""
|
||||
|
||||
def __init__(self, channel_energies: Optional[np.ndarray] = None):
|
||||
self._channel_energies = channel_energies
|
||||
self._calibration: Optional[CalibrationResult] = None
|
||||
self._calibration_points: Optional[list[tuple[float, float]]] = None
|
||||
|
||||
@property
|
||||
def is_calibrated(self) -> bool:
|
||||
"""Whether a calibration has been applied."""
|
||||
return self._calibration is not None
|
||||
|
||||
def add_calibration_point(
|
||||
self,
|
||||
channel: int | float,
|
||||
energy_keV: float,
|
||||
) -> None:
|
||||
"""Add a calibration point (known line energy → channel number).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
channel : int or float
|
||||
Detector channel number.
|
||||
energy_keV : float
|
||||
Known energy of the calibration line (keV).
|
||||
"""
|
||||
if self._calibration_points is None:
|
||||
self._calibration_points = []
|
||||
self._calibration_points.append((float(channel), float(energy_keV)))
|
||||
|
||||
def set_calibration_points(self, points: list[tuple[float, float]]) -> None:
|
||||
"""Set calibration points directly.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
points : list of (channel, energy_keV)
|
||||
Known calibration line positions.
|
||||
"""
|
||||
self._calibration_points = [(float(c), float(e)) for c, e in points]
|
||||
|
||||
def fit_calibration(
|
||||
self,
|
||||
order: int = 1,
|
||||
channel_energies: Optional[np.ndarray] = None,
|
||||
) -> CalibrationResult:
|
||||
"""Fit calibration polynomial using known points.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
order : int
|
||||
Polynomial order (1 = linear, 2 = quadratic).
|
||||
channel_energies : np.ndarray, optional
|
||||
Full channel-to-energy mapping. If provided, fits to
|
||||
all points; otherwise fits only to added calibration points.
|
||||
|
||||
Returns
|
||||
-------
|
||||
CalibrationResult
|
||||
"""
|
||||
if channel_energies is not None:
|
||||
# Fit to full mapping
|
||||
channels = np.arange(len(channel_energies), dtype=np.float64)
|
||||
x = channels
|
||||
y = np.asarray(channel_energies, dtype=np.float64)
|
||||
elif self._calibration_points is not None:
|
||||
x = np.array([p[0] for p in self._calibration_points])
|
||||
y = np.array([p[1] for p in self._calibration_points])
|
||||
else:
|
||||
raise ValueError(
|
||||
"No calibration data available. Use add_calibration_point() "
|
||||
"or pass channel_energies."
|
||||
)
|
||||
|
||||
# Fit polynomial: E(ch) = a₀ + a₁·ch + a₂·ch² + ...
|
||||
try:
|
||||
coeffs, cov = np.polyfit(x, y, order, cov=True)
|
||||
except ValueError:
|
||||
# Too few points for covariance estimation
|
||||
coeffs = np.polyfit(x, y, order)
|
||||
cov = np.eye(order + 1) * 1e-10
|
||||
|
||||
# Compute χ²
|
||||
y_fit = np.polyval(coeffs[::-1], x)
|
||||
# Assume uniform uncertainty of 0.1 keV for calibration points
|
||||
uncertainties = np.full_like(x, 0.1)
|
||||
chi2 = float(np.sum(((y - y_fit) / uncertainties) ** 2))
|
||||
n_dof = max(len(x) - order - 1, 1)
|
||||
|
||||
return CalibrationResult(
|
||||
coefficients=coeffs,
|
||||
order=order,
|
||||
chi2_per_dof=chi2 / n_dof,
|
||||
n_points=len(x),
|
||||
covariance=cov,
|
||||
)
|
||||
|
||||
def apply_calibration(
|
||||
self,
|
||||
order: int = 1,
|
||||
channel_energies: Optional[np.ndarray] = None,
|
||||
) -> CalibrationResult:
|
||||
"""Fit and store calibration for future use.
|
||||
|
||||
See :meth:`fit_calibration` for parameters.
|
||||
"""
|
||||
self._calibration = self.fit_calibration(order, channel_energies)
|
||||
return self._calibration
|
||||
|
||||
def channel_to_energy(self, channel: float | np.ndarray) -> np.ndarray:
|
||||
"""Convert channel number to energy using current calibration.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
channel : float or np.ndarray
|
||||
Channel number(s).
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Energy in keV.
|
||||
|
||||
Raises
|
||||
------
|
||||
RuntimeError
|
||||
If no calibration has been applied.
|
||||
"""
|
||||
if self._calibration is None:
|
||||
raise RuntimeError(
|
||||
"No calibration applied. Call apply_calibration() first."
|
||||
)
|
||||
return self._calibration(channel)
|
||||
|
||||
def energy_to_channel(self, energy_keV: float | np.ndarray) -> np.ndarray:
|
||||
"""Convert energy to channel using current calibration.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
energy_keV : float or np.ndarray
|
||||
Energy in keV.
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Channel number(s).
|
||||
|
||||
Raises
|
||||
------
|
||||
RuntimeError
|
||||
If no calibration has been applied.
|
||||
"""
|
||||
if self._calibration is None:
|
||||
raise RuntimeError(
|
||||
"No calibration applied. Call apply_calibration() first."
|
||||
)
|
||||
inv = self._calibration.invert()
|
||||
return inv(energy_keV)
|
||||
|
||||
def get_linear_calibration(
|
||||
self,
|
||||
low_energy_keV: float,
|
||||
high_energy_keV: float,
|
||||
low_channel: int,
|
||||
high_channel: int,
|
||||
) -> CalibrationResult:
|
||||
"""Create a linear calibration from two known points.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
low_energy_keV : float
|
||||
Energy of first calibration line (keV).
|
||||
high_energy_keV : float
|
||||
Energy of second calibration line (keV).
|
||||
low_channel : int
|
||||
Channel number of first line.
|
||||
high_channel : int
|
||||
Channel number of second line.
|
||||
|
||||
Returns
|
||||
-------
|
||||
CalibrationResult
|
||||
"""
|
||||
self.set_calibration_points(
|
||||
[
|
||||
(float(low_channel), low_energy_keV),
|
||||
(float(high_channel), high_energy_keV),
|
||||
]
|
||||
)
|
||||
return self.apply_calibration(order=1)
|
||||
@@ -0,0 +1,268 @@
|
||||
"""
|
||||
Instrumental corrections for experimental spectra.
|
||||
|
||||
Provides dead-time correction, pile-up correction, and background
|
||||
subtraction for beta spectrometry data.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from typing import Optional
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
class DeadTimeCorrection:
|
||||
"""Dead-time correction for paralyzable or non-paralyzable detectors.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
dead_time : float
|
||||
Dead time per event in seconds.
|
||||
model : str
|
||||
'non_paralyzable' (default) or 'paralyzable'.
|
||||
"""
|
||||
|
||||
def __init__(self, dead_time: float, model: str = "non_paralyzable"):
|
||||
self.dead_time = dead_time
|
||||
self.model = model
|
||||
assert model in (
|
||||
"non_paralyzable",
|
||||
"paralyzable",
|
||||
), f"Unknown dead-time model: {model}"
|
||||
assert dead_time >= 0, "Dead time must be non-negative"
|
||||
|
||||
def live_time_fraction(self, observed_rate: float) -> float:
|
||||
"""Calculate live time fraction given observed count rate.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
observed_rate : float
|
||||
Observed count rate in counts/second.
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Fraction of time the detector is live (0 to 1).
|
||||
"""
|
||||
if self.model == "non_paralyzable":
|
||||
return max(0.0, 1.0 - observed_rate * self.dead_time)
|
||||
else: # paralyzable
|
||||
return np.exp(-observed_rate * self.dead_time)
|
||||
|
||||
def correct_rate(self, observed_rate: float) -> float:
|
||||
"""Correct observed rate for dead time.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
observed_rate : float
|
||||
Observed count rate (counts/s).
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
True count rate (counts/s).
|
||||
"""
|
||||
live_frac = self.live_time_fraction(observed_rate)
|
||||
if live_frac <= 0:
|
||||
return float("inf")
|
||||
return observed_rate / live_frac
|
||||
|
||||
def correct_counts(self, counts: np.ndarray, live_time: float) -> np.ndarray:
|
||||
"""Correct counts for dead time.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
counts : np.ndarray
|
||||
Measured counts per bin.
|
||||
live_time : float
|
||||
Total live time in seconds.
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Dead-time corrected counts.
|
||||
"""
|
||||
total_observed = np.sum(counts)
|
||||
total_true = self.correct_rate(total_observed / live_time) * live_time
|
||||
if total_observed <= 0:
|
||||
return counts.copy()
|
||||
return counts * (total_true / total_observed)
|
||||
|
||||
|
||||
class PileUpCorrection:
|
||||
"""Simple pile-up correction for beta spectrometry.
|
||||
|
||||
Estimates the fraction of pile-up events based on observed count rate
|
||||
and dead time, and redistributes pile-up events back into the spectrum.
|
||||
|
||||
This is a first-order approximation. For precise pile-up correction,
|
||||
digital pulse processing or anti-pile-up electronics are recommended.
|
||||
"""
|
||||
|
||||
def __init__(self, dead_time: float):
|
||||
self.dead_time = dead_time
|
||||
|
||||
def estimate_pile_up_fraction(self, count_rate: float) -> float:
|
||||
"""Estimate pile-up fraction from count rate.
|
||||
|
||||
For a non-paralyzable detector, pile-up fraction ≈ (count_rate × dead_time)².
|
||||
|
||||
Parameters
|
||||
----------
|
||||
count_rate : float
|
||||
Observed count rate in counts/second.
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Estimated fraction of events affected by pile-up.
|
||||
"""
|
||||
return (count_rate * self.dead_time) ** 2
|
||||
|
||||
def correct_spectrum(
|
||||
self,
|
||||
counts: np.ndarray,
|
||||
live_time: float,
|
||||
background: Optional[np.ndarray] = None,
|
||||
) -> np.ndarray:
|
||||
"""Apply first-order pile-up correction.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
counts : np.ndarray
|
||||
Measured counts per bin.
|
||||
live_time : float
|
||||
Live time in seconds.
|
||||
background : np.ndarray, optional
|
||||
Background counts per bin.
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Pile-up corrected counts.
|
||||
"""
|
||||
if background is not None:
|
||||
net_counts = counts - background
|
||||
else:
|
||||
net_counts = counts.copy()
|
||||
|
||||
observed_rate = np.sum(net_counts) / live_time
|
||||
pileup_frac = self.estimate_pile_up_fraction(observed_rate)
|
||||
|
||||
# Redistribute pile-up events proportionally to the spectrum shape
|
||||
if pileup_frac > 0 and np.sum(net_counts) > 0:
|
||||
correction = net_counts * pileup_frac / (1.0 - pileup_frac)
|
||||
return net_counts + correction
|
||||
|
||||
return net_counts
|
||||
|
||||
|
||||
class BackgroundSubtractor:
|
||||
"""Background subtraction for beta spectra.
|
||||
|
||||
Supports:
|
||||
- Constant background (from sidebands)
|
||||
- Polynomial background
|
||||
- Region-of-interest based background estimation
|
||||
"""
|
||||
|
||||
@staticmethod
|
||||
def constant_background(
|
||||
counts: np.ndarray,
|
||||
energies: np.ndarray,
|
||||
low_roi: tuple[float, float],
|
||||
high_roi: tuple[float, float],
|
||||
) -> np.ndarray:
|
||||
"""Subtract constant background estimated from sidebands.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
counts : np.ndarray
|
||||
Measured counts per bin.
|
||||
energies : np.ndarray
|
||||
Energy values (keV) per bin.
|
||||
low_roi : (E_min, E_max)
|
||||
Energy range for low-side background estimation.
|
||||
high_roi : (E_min, E_max)
|
||||
Energy range for high-side background estimation.
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Counts with constant background subtracted.
|
||||
"""
|
||||
low_mask = (energies >= low_roi[0]) & (energies <= low_roi[1])
|
||||
high_mask = (energies >= high_roi[0]) & (energies <= high_roi[1])
|
||||
|
||||
low_bg = np.mean(counts[low_mask]) if np.any(low_mask) else 0.0
|
||||
high_bg = np.mean(counts[high_mask]) if np.any(high_mask) else 0.0
|
||||
constant_bg = (low_bg + high_bg) / 2.0
|
||||
|
||||
return counts - constant_bg
|
||||
|
||||
@staticmethod
|
||||
def polynomial_background(
|
||||
counts: np.ndarray,
|
||||
energies: np.ndarray,
|
||||
order: int = 2,
|
||||
low_roi: tuple[float, float] = (0.0, 0.0),
|
||||
high_roi: tuple[float, float] = (0.0, 0.0),
|
||||
) -> np.ndarray:
|
||||
"""Subtract polynomial background fitted from sidebands.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
counts : np.ndarray
|
||||
Measured counts per bin.
|
||||
energies : np.ndarray
|
||||
Energy values (keV) per bin.
|
||||
order : int
|
||||
Polynomial order for background fit.
|
||||
low_roi : (E_min, E_max)
|
||||
Energy range for low-side background fitting.
|
||||
high_roi : (E_min, E_max)
|
||||
Energy range for high-side background fitting.
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Counts with polynomial background subtracted.
|
||||
"""
|
||||
low_mask = (energies >= low_roi[0]) & (energies <= low_roi[1])
|
||||
high_mask = (energies >= high_roi[0]) & (energies <= high_roi[1])
|
||||
|
||||
mask = low_mask | high_mask
|
||||
if not np.any(mask):
|
||||
return counts
|
||||
|
||||
x_fit = energies[mask]
|
||||
y_fit = counts[mask]
|
||||
|
||||
coeffs = np.polyfit(x_fit, y_fit, order)
|
||||
bg = np.polyval(coeffs[::-1], energies)
|
||||
|
||||
return counts - bg
|
||||
|
||||
@staticmethod
|
||||
def roi_average(
|
||||
counts: np.ndarray,
|
||||
roi_mask: np.ndarray,
|
||||
) -> float:
|
||||
"""Calculate average background from a region of interest.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
counts : np.ndarray
|
||||
Measured counts per bin.
|
||||
roi_mask : np.ndarray
|
||||
Boolean mask selecting background region.
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Average background count per bin.
|
||||
"""
|
||||
if not np.any(roi_mask):
|
||||
return 0.0
|
||||
return float(np.mean(counts[roi_mask]))
|
||||
@@ -0,0 +1,307 @@
|
||||
"""
|
||||
Simple peak fitters for calibration source analysis.
|
||||
|
||||
Provides Gaussian and Gaussian+background fitting for identifying
|
||||
known X-ray and gamma-ray lines used in energy calibration.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from typing import Optional, Tuple
|
||||
|
||||
import numpy as np
|
||||
from scipy.optimize import least_squares
|
||||
|
||||
|
||||
@dataclass
|
||||
class PeakFitResult:
|
||||
"""Result of a peak fit."""
|
||||
|
||||
center_keV: float
|
||||
"""Peak center energy in keV."""
|
||||
|
||||
sigma_keV: float
|
||||
"""Peak width (σ) in keV."""
|
||||
|
||||
amplitude: float
|
||||
"""Peak amplitude (counts at peak)."""
|
||||
|
||||
background: float
|
||||
"""Constant background level."""
|
||||
|
||||
chi2: float
|
||||
"""Reduced χ²."""
|
||||
|
||||
n_dof: int
|
||||
"""Number of degrees of freedom."""
|
||||
|
||||
covariance: np.ndarray
|
||||
"""Parameter covariance matrix.
|
||||
|
||||
Parameter order: [amplitude, center, sigma, background]
|
||||
"""
|
||||
|
||||
success: bool
|
||||
"""Whether the fit converged."""
|
||||
|
||||
def resolution(self) -> float:
|
||||
"""FWHM in keV."""
|
||||
return 2.35482 * self.sigma_keV
|
||||
|
||||
def resolution_fraction(self) -> float:
|
||||
"""FWHM / center as a fraction."""
|
||||
if self.center_keV > 0:
|
||||
return self.resolution() / self.center_keV
|
||||
return 0.0
|
||||
|
||||
|
||||
class GaussianFitter:
|
||||
"""Fit a Gaussian peak on top of constant background.
|
||||
|
||||
Model: f(x) = A · exp(-(x - μ)² / (2σ²)) + B
|
||||
|
||||
Parameters
|
||||
----------
|
||||
energies : np.ndarray
|
||||
Energy values (keV), shape (N,).
|
||||
counts : np.ndarray
|
||||
Measured counts, shape (N,).
|
||||
uncertainties : np.ndarray, optional
|
||||
Measurement uncertainties. If None, assumed to be 1.0.
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
energies: np.ndarray,
|
||||
counts: np.ndarray,
|
||||
uncertainties: Optional[np.ndarray] = None,
|
||||
):
|
||||
self.energies = np.asarray(energies, dtype=np.float64)
|
||||
self.counts = np.asarray(counts, dtype=np.float64)
|
||||
|
||||
if uncertainties is not None:
|
||||
self.uncertainties = np.asarray(uncertainties, dtype=np.float64)
|
||||
else:
|
||||
self.uncertainties = np.ones_like(self.counts)
|
||||
|
||||
# Filter out zero-uncertainty points
|
||||
valid = self.uncertainties > 0
|
||||
self.energies = self.energies[valid]
|
||||
self.counts = self.counts[valid]
|
||||
self.uncertainties = self.uncertainties[valid]
|
||||
|
||||
def _gaussian_background(self, params: np.ndarray, x: np.ndarray) -> np.ndarray:
|
||||
"""Gaussian + constant background model."""
|
||||
A, mu, sigma, B = params
|
||||
return A * np.exp(-0.5 * ((x - mu) / sigma) ** 2) + B
|
||||
|
||||
def _residuals(self, params: np.ndarray) -> np.ndarray:
|
||||
"""Weighted residuals."""
|
||||
model = self._gaussian_background(params, self.energies)
|
||||
return (self.counts - model) / self.uncertainties
|
||||
|
||||
def _chi2(self, params: np.ndarray) -> float:
|
||||
"""Compute χ²."""
|
||||
residuals = self._residuals(params)
|
||||
return float(np.sum(residuals**2))
|
||||
|
||||
def _compute_covariance(self, params: np.ndarray, jac: np.ndarray) -> np.ndarray:
|
||||
"""Estimate covariance from Jacobian."""
|
||||
jtj = jac.T @ jac
|
||||
if np.linalg.cond(jtj) < 1e12:
|
||||
cov = np.linalg.inv(jtj)
|
||||
chi2 = self._chi2(params)
|
||||
ndof = max(len(self.energies) - 4, 1)
|
||||
cov = cov * chi2 / ndof
|
||||
return cov
|
||||
return np.diag([1e10, 1e6, 1e4, 1e6])
|
||||
|
||||
def fit(
|
||||
self,
|
||||
guess: Optional[Tuple[float, float, float, float]] = None,
|
||||
bounds: Optional[Tuple[tuple, tuple]] = None,
|
||||
) -> PeakFitResult:
|
||||
"""Fit a Gaussian peak to the data.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
guess : (amplitude, center, sigma, background), optional
|
||||
Initial parameter guess. Auto-estimated if None.
|
||||
bounds : ((A_min, μ_min, σ_min, B_min), (A_max, μ_max, σ_max, B_max)), optional
|
||||
Parameter bounds.
|
||||
|
||||
Returns
|
||||
-------
|
||||
PeakFitResult
|
||||
"""
|
||||
if guess is None:
|
||||
guess = self._estimate_guess()
|
||||
|
||||
if bounds is None:
|
||||
bounds = (
|
||||
(0.0, min(self.energies), 0.1, 0.0),
|
||||
(np.inf, max(self.energies), 50.0, np.inf),
|
||||
)
|
||||
|
||||
result = least_squares(
|
||||
self._residuals,
|
||||
x0=guess,
|
||||
bounds=bounds,
|
||||
method="trf",
|
||||
max_nfev=1000,
|
||||
)
|
||||
|
||||
A, mu, sigma, B = result.x
|
||||
n_dof = max(len(self.energies) - 4, 1)
|
||||
chi2 = self._chi2(result.x)
|
||||
|
||||
# Compute covariance from Jacobian
|
||||
_, jac = result.fun, result.jac
|
||||
cov = self._compute_covariance(result.x, jac) if jac is not None else np.eye(4)
|
||||
|
||||
return PeakFitResult(
|
||||
center_keV=mu,
|
||||
sigma_keV=sigma,
|
||||
amplitude=A,
|
||||
background=B,
|
||||
chi2=chi2 / n_dof if n_dof > 0 else chi2,
|
||||
n_dof=n_dof,
|
||||
covariance=cov,
|
||||
success=result.success,
|
||||
)
|
||||
|
||||
def _estimate_guess(self) -> Tuple[float, float, float, float]:
|
||||
"""Auto-estimate initial parameters from data."""
|
||||
# Find the peak
|
||||
peak_idx = np.argmax(self.counts)
|
||||
center = self.energies[peak_idx]
|
||||
amplitude = self.counts[peak_idx]
|
||||
|
||||
# Estimate sigma from FWHM
|
||||
half_max = amplitude / 2.0
|
||||
left_idx = np.argmax(self.counts[:peak_idx] < half_max)
|
||||
right_idx = np.argmin(self.counts[peak_idx:] < half_max) + peak_idx
|
||||
|
||||
if left_idx < peak_idx and right_idx > peak_idx:
|
||||
fwhm = self.energies[right_idx] - self.energies[left_idx]
|
||||
sigma = fwhm / 2.35482
|
||||
else:
|
||||
sigma = 1.0 # Default width
|
||||
|
||||
# Background estimate from average of edges
|
||||
n_edge = max(10, len(self.counts) // 20)
|
||||
background = np.mean(
|
||||
np.concatenate([self.counts[:n_edge], self.counts[-n_edge:]])
|
||||
)
|
||||
|
||||
return (amplitude - background, center, max(sigma, 0.1), background)
|
||||
|
||||
|
||||
class PeakFitter:
|
||||
"""Find and fit multiple peaks in a spectrum.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
energies : np.ndarray
|
||||
Energy values (keV).
|
||||
counts : np.ndarray
|
||||
Measured counts.
|
||||
uncertainties : np.ndarray, optional
|
||||
Measurement uncertainties.
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
energies: np.ndarray,
|
||||
counts: np.ndarray,
|
||||
uncertainties: Optional[np.ndarray] = None,
|
||||
):
|
||||
self.energies = np.asarray(energies, dtype=np.float64)
|
||||
self.counts = np.asarray(counts, dtype=np.float64)
|
||||
|
||||
if uncertainties is not None:
|
||||
self.uncertainties = np.asarray(uncertainties, dtype=np.float64)
|
||||
else:
|
||||
self.uncertainties = np.sqrt(np.maximum(self.counts, 1.0))
|
||||
|
||||
def find_peaks(
|
||||
self,
|
||||
min_height: float = 10.0,
|
||||
min_distance_keV: float = 2.0,
|
||||
smoothing_window: int = 5,
|
||||
) -> np.ndarray:
|
||||
"""Find candidate peak positions.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
min_height : float
|
||||
Minimum peak height above background (counts).
|
||||
min_distance_keV : float
|
||||
Minimum separation between peaks in keV.
|
||||
smoothing_window : int
|
||||
Smoothing window for peak detection (odd number).
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Energy values of detected peak centers.
|
||||
"""
|
||||
from scipy.signal import find_peaks
|
||||
|
||||
# Estimate background by smoothing
|
||||
window = smoothing_window if smoothing_window % 2 == 1 else smoothing_window + 1
|
||||
if len(self.counts) > window:
|
||||
from scipy.ndimage import uniform_filter1d
|
||||
|
||||
background_est = uniform_filter1d(self.counts, size=window)
|
||||
else:
|
||||
background_est = np.full_like(self.counts, np.mean(self.counts))
|
||||
|
||||
heights = self.counts - background_est
|
||||
indices, _ = find_peaks(
|
||||
heights,
|
||||
height=min_height,
|
||||
distance=int(min_distance_keV / max(1.0, np.diff(self.energies).min())),
|
||||
)
|
||||
|
||||
return self.energies[indices]
|
||||
|
||||
def fit_peaks(
|
||||
self,
|
||||
peak_energies: Optional[np.ndarray] = None,
|
||||
known_energies: Optional[np.ndarray] = None,
|
||||
) -> list[PeakFitResult]:
|
||||
"""Fit peaks at given or auto-detected positions.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
peak_energies : np.ndarray, optional
|
||||
Explicit peak center energies to fit.
|
||||
known_energies : np.ndarray, optional
|
||||
Known calibration line energies — auto-detect peaks near these.
|
||||
|
||||
Returns
|
||||
-------
|
||||
list[PeakFitResult]
|
||||
"""
|
||||
if known_energies is not None:
|
||||
# Find peaks near known energies
|
||||
peak_energies = np.array([self._nearest_channel(e) for e in known_energies])
|
||||
|
||||
if peak_energies is None:
|
||||
peak_energies = self.find_peaks()
|
||||
|
||||
results = []
|
||||
for center in peak_energies:
|
||||
fitter = GaussianFitter(self.energies, self.counts, self.uncertainties)
|
||||
result = fitter.fit(guess=(10.0, center, 1.0, 0.0))
|
||||
if result.success:
|
||||
results.append(result)
|
||||
|
||||
return results
|
||||
|
||||
def _nearest_channel(self, target_keV: float) -> float:
|
||||
"""Find the energy channel closest to target_keV."""
|
||||
idx = np.argmin(np.abs(self.energies - target_keV))
|
||||
return self.energies[idx]
|
||||
@@ -0,0 +1,161 @@
|
||||
"""
|
||||
Experimental spectrum container.
|
||||
|
||||
Stores measured detector spectra with metadata, uncertainties,
|
||||
and calibration information.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
from typing import Optional
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
@dataclass
|
||||
class ExpSpectrum:
|
||||
"""Container for an experimental beta spectrum measurement.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
energies : np.ndarray
|
||||
Energy calibration values for each channel (keV), shape (N_channels,).
|
||||
counts : np.ndarray
|
||||
Raw counts per channel, shape (N_channels,).
|
||||
errors : np.ndarray, optional
|
||||
Statistical uncertainties per channel. Defaults to sqrt(counts).
|
||||
metadata : dict, optional
|
||||
Arbitrary metadata (source, run ID, date, etc.).
|
||||
dead_time : float, optional
|
||||
Dead time fraction (0.0 to 1.0).
|
||||
live_time : float, optional
|
||||
Live acquisition time in seconds.
|
||||
source : str, optional
|
||||
Radioactive source identifier.
|
||||
run_id : str, optional
|
||||
Unique run identifier.
|
||||
date : str, optional
|
||||
Measurement date (ISO format).
|
||||
"""
|
||||
|
||||
energies: np.ndarray
|
||||
counts: np.ndarray
|
||||
errors: Optional[np.ndarray] = None
|
||||
metadata: dict = field(default_factory=dict)
|
||||
dead_time: float = 0.0
|
||||
live_time: Optional[float] = None
|
||||
source: Optional[str] = None
|
||||
run_id: Optional[str] = None
|
||||
date: Optional[str] = None
|
||||
|
||||
def __post_init__(self) -> None:
|
||||
"""Validate and compute default errors."""
|
||||
self.energies = np.asarray(self.energies, dtype=np.float64)
|
||||
self.counts = np.asarray(self.counts, dtype=np.float64)
|
||||
assert self.energies.shape == self.counts.shape, (
|
||||
f"Energies ({self.energies.shape}) and counts ({self.counts.shape}) "
|
||||
"must have the same shape"
|
||||
)
|
||||
assert np.all(self.energies >= 0), "Energies must be non-negative"
|
||||
assert np.all(self.counts >= 0), "Counts must be non-negative"
|
||||
|
||||
if self.errors is None:
|
||||
# Poisson statistics: σ = √N (handle zero counts gracefully)
|
||||
self.errors = np.sqrt(np.maximum(self.counts, 1.0))
|
||||
else:
|
||||
self.errors = np.asarray(self.errors, dtype=np.float64)
|
||||
assert self.errors.shape == self.counts.shape
|
||||
assert np.all(self.errors > 0), "Errors must be positive"
|
||||
|
||||
@property
|
||||
def n_channels(self) -> int:
|
||||
"""Number of energy channels."""
|
||||
return len(self.energies)
|
||||
|
||||
@property
|
||||
def energy_range(self) -> tuple[float, float]:
|
||||
"""(min_energy_keV, max_energy_keV)."""
|
||||
return (float(self.energies[0]), float(self.energies[-1]))
|
||||
|
||||
@property
|
||||
def total_counts(self) -> float:
|
||||
"""Total counts in the spectrum."""
|
||||
return float(np.sum(self.counts))
|
||||
|
||||
@property
|
||||
def live_time_corrected_rate(self) -> Optional[float]:
|
||||
"""Count rate corrected for dead time, in counts per second."""
|
||||
if self.live_time is None or self.live_time <= 0:
|
||||
return None
|
||||
live_seconds = self.live_time * (1.0 - self.dead_time)
|
||||
return self.total_counts / live_seconds
|
||||
|
||||
def normalize(self) -> np.ndarray:
|
||||
"""Return spectrum normalized to unit area."""
|
||||
integral = np.trapezoid(self.counts, self.energies)
|
||||
if integral > 0:
|
||||
return self.counts / integral
|
||||
return self.counts.copy()
|
||||
|
||||
def apply_dead_time_correction(self) -> ExpSpectrum:
|
||||
"""Return a new spectrum corrected for dead time.
|
||||
|
||||
Correction: N_true = N_measured / (1 - τ), where τ is dead time fraction.
|
||||
"""
|
||||
if self.dead_time <= 0 or self.dead_time >= 1.0:
|
||||
return ExpSpectrum(
|
||||
energies=self.energies.copy(),
|
||||
counts=self.counts.copy(),
|
||||
errors=self.errors.copy(),
|
||||
metadata=self.metadata.copy(),
|
||||
dead_time=self.dead_time,
|
||||
live_time=self.live_time,
|
||||
source=self.source,
|
||||
run_id=self.run_id,
|
||||
date=self.date,
|
||||
)
|
||||
|
||||
corrected_counts = self.counts / (1.0 - self.dead_time)
|
||||
corrected_errors = self.errors / (1.0 - self.dead_time)
|
||||
|
||||
return ExpSpectrum(
|
||||
energies=self.energies.copy(),
|
||||
counts=corrected_counts,
|
||||
errors=corrected_errors,
|
||||
metadata=self.metadata.copy(),
|
||||
dead_time=0.0, # Corrected
|
||||
live_time=self.live_time,
|
||||
source=self.source,
|
||||
run_id=self.run_id,
|
||||
date=self.date,
|
||||
)
|
||||
|
||||
def to_dict(self) -> dict:
|
||||
"""Serialize to dictionary."""
|
||||
return {
|
||||
"energies": self.energies.tolist(),
|
||||
"counts": self.counts.tolist(),
|
||||
"errors": self.errors.tolist(),
|
||||
"metadata": self.metadata,
|
||||
"dead_time": self.dead_time,
|
||||
"live_time": self.live_time,
|
||||
"source": self.source,
|
||||
"run_id": self.run_id,
|
||||
"date": self.date,
|
||||
}
|
||||
|
||||
@classmethod
|
||||
def from_dict(cls, data: dict) -> ExpSpectrum:
|
||||
"""Deserialize from dictionary."""
|
||||
return cls(
|
||||
energies=np.array(data["energies"], dtype=np.float64),
|
||||
counts=np.array(data["counts"], dtype=np.float64),
|
||||
errors=np.array(data["errors"], dtype=np.float64),
|
||||
metadata=data.get("metadata", {}),
|
||||
dead_time=data.get("dead_time", 0.0),
|
||||
live_time=data.get("live_time"),
|
||||
source=data.get("source"),
|
||||
run_id=data.get("run_id"),
|
||||
date=data.get("date"),
|
||||
)
|
||||
@@ -0,0 +1,38 @@
|
||||
"""
|
||||
Fitter module — main analysis framework for beta spectrometry.
|
||||
|
||||
Provides the analysis workhorse: theoretical model + detector convolution,
|
||||
multi-parameter fitting, and physics parameter extraction (C(W), g_V, g_A).
|
||||
|
||||
Usage:
|
||||
from fitter import SpectrumModel, SpectrumFitter, AnalysisResult
|
||||
|
||||
# Build theoretical model with detector response
|
||||
model = SpectrumModel(spectrum, detector_response, channel_energies)
|
||||
|
||||
# Fit model to experimental data
|
||||
fitter = SpectrumFitter(model, exp_spectrum)
|
||||
result = fitter.fit(
|
||||
x0=[endpoint_keV, normalization, background],
|
||||
bounds=[(...), (...), (...)],
|
||||
)
|
||||
|
||||
# Extract physics parameters
|
||||
print(result.summary())
|
||||
"""
|
||||
|
||||
from .model import SpectrumModel
|
||||
from .fit_engine import SpectrumFitter, CurveFitter, FitConfig
|
||||
from .result import AnalysisResult, FitSummary
|
||||
from .extractor import CWExtractor, GVAExtractor
|
||||
|
||||
__all__ = [
|
||||
"SpectrumModel",
|
||||
"SpectrumFitter",
|
||||
"CurveFitter",
|
||||
"FitConfig",
|
||||
"AnalysisResult",
|
||||
"FitSummary",
|
||||
"CWExtractor",
|
||||
"GVAExtractor",
|
||||
]
|
||||
@@ -0,0 +1,373 @@
|
||||
"""
|
||||
C(W) shape factor and g_V/g_A extraction from fitted spectra.
|
||||
|
||||
Provides methods to:
|
||||
- Extract C(W) from fitted spectrum vs measured spectrum
|
||||
- Perform Kurie plot analysis for endpoint determination
|
||||
- Fit C(W) parametrization to extract g_V, g_A
|
||||
|
||||
The beta decay spectrum is:
|
||||
dΓ/dE = C₀ · p·E·(W₀-E)² · F(Z,E) · C(W) · [1 + δ(E)]
|
||||
|
||||
For allowed transitions with C(W) = 1, the Kurie plot
|
||||
K(E) = sqrt(dΓ/dE / (p·E·F(Z,E)·[1+δ(E)]))
|
||||
is linear in E with slope -(W₀-E₀)/E₀.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
from typing import TYPE_CHECKING, Any, Optional, Tuple
|
||||
|
||||
import numpy as np
|
||||
|
||||
if TYPE_CHECKING:
|
||||
from beta_spectrum.spectrum import SpectrumConfig
|
||||
|
||||
from beta_spectrum.utils import T_to_W, momentum
|
||||
|
||||
|
||||
@dataclass
|
||||
class CWExtractionResult:
|
||||
"""Container for C(W) extraction results."""
|
||||
|
||||
energies_W: np.ndarray
|
||||
"""Total energy grid in m_e units."""
|
||||
|
||||
energies_keV: np.ndarray
|
||||
"""Kinetic energy grid in keV."""
|
||||
|
||||
cw_values: np.ndarray
|
||||
"""Extracted C(W) values."""
|
||||
|
||||
cw_errors: np.ndarray
|
||||
"""Uncertainties on C(W) values."""
|
||||
|
||||
flux: np.ndarray
|
||||
"""Measured flux (counts per energy bin)."""
|
||||
|
||||
flux_errors: np.ndarray
|
||||
"""Uncertainties on flux."""
|
||||
|
||||
endpoint_W: float
|
||||
"""Endpoint energy W₀ in m_e units."""
|
||||
|
||||
endpoint_keV: float
|
||||
"""Endpoint kinetic energy in keV."""
|
||||
|
||||
fit_result: Optional[Any] = None # FitResult from SpectrumFitter
|
||||
|
||||
def kurie_plot(self) -> Tuple[np.ndarray, np.ndarray]:
|
||||
"""Compute Kurie plot values.
|
||||
|
||||
Returns
|
||||
-------
|
||||
(energies_keV, kurie_values)
|
||||
"""
|
||||
p = momentum(self.energies_W)
|
||||
E = self.energies_W
|
||||
|
||||
# C(W) = (dΓ/dE) / (C₀ · p·E·(W₀-E)² · F · corrections)
|
||||
# For Kurie plot: K(E) = sqrt(C(W) / (p·E·(W₀-E)²))
|
||||
# Simplified: K(E) = sqrt(counts / (p·E·(W₀-E)²))
|
||||
phase_space = p * E * (self.endpoint_W - E) ** 2
|
||||
mask = phase_space > 0
|
||||
|
||||
kurie = np.zeros_like(self.energies_keV)
|
||||
kurie[mask] = np.sqrt(
|
||||
self.flux[mask] / (phase_space[mask] * self.energies_keV[mask] ** 2)
|
||||
)
|
||||
|
||||
return self.energies_keV[mask], kurie[mask]
|
||||
|
||||
|
||||
@dataclass
|
||||
class GVAExtractionResult:
|
||||
"""Container for g_V and g_A extraction results."""
|
||||
|
||||
gv_eff: float
|
||||
"""Effective vector coupling constant."""
|
||||
|
||||
gv_error: float
|
||||
"""1σ uncertainty on g_V."""
|
||||
|
||||
ga_eff: float
|
||||
"""Effective axial coupling constant."""
|
||||
|
||||
ga_error: float
|
||||
"""1σ uncertainty on g_A."""
|
||||
|
||||
correlation: float
|
||||
"""Correlation coefficient between g_V and g_A."""
|
||||
|
||||
chi2_per_dof: float
|
||||
"""Reduced χ²."""
|
||||
|
||||
def summary(self) -> str:
|
||||
lines = [
|
||||
"=" * 60,
|
||||
"g_V / g_A Extraction",
|
||||
"=" * 60,
|
||||
f"g_V^eff = {self.gv_eff:.6f} ± {self.gv_error:.6f}",
|
||||
f"g_A^eff = {self.ga_eff:.6f} ± {self.ga_error:.6f}",
|
||||
f"Correlation(g_V, g_A) = {self.correlation:.3f}",
|
||||
f"χ²/ndof = {self.chi2_per_dof:.3f}",
|
||||
"=" * 60,
|
||||
]
|
||||
return "\n".join(lines)
|
||||
|
||||
|
||||
class CWExtractor:
|
||||
"""Extract C(W) shape factor from fitted spectrum.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
config : SpectrumConfig
|
||||
Configuration used for the theoretical model.
|
||||
"""
|
||||
|
||||
def __init__(self, config: "SpectrumConfig"):
|
||||
self.config = config
|
||||
|
||||
def extract(
|
||||
self,
|
||||
measured_counts: np.ndarray,
|
||||
measured_errors: np.ndarray,
|
||||
energies_keV: np.ndarray,
|
||||
model_counts: np.ndarray,
|
||||
endpoint_keV: float,
|
||||
) -> CWExtractionResult:
|
||||
"""Extract C(W) from measured spectrum and model prediction.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
measured_counts : np.ndarray
|
||||
Measured counts per channel.
|
||||
measured_errors : np.ndarray
|
||||
Measurement uncertainties.
|
||||
energies_keV : np.ndarray
|
||||
Energy grid in keV.
|
||||
model_counts : np.ndarray
|
||||
Theoretical model counts (from fit).
|
||||
endpoint_keV : float
|
||||
Endpoint energy in keV.
|
||||
|
||||
Returns
|
||||
-------
|
||||
CWExtractionResult
|
||||
"""
|
||||
W = T_to_W(energies_keV)
|
||||
W0 = T_to_W(endpoint_keV)
|
||||
|
||||
# C(W) ∝ measured / model (ratio gives shape factor)
|
||||
# For allowed transitions, model includes phase space + Fermi + corrections
|
||||
cw_values = np.where(model_counts > 0, measured_counts / model_counts, 1.0)
|
||||
cw_errors = np.where(
|
||||
model_counts > 0,
|
||||
measured_errors / model_counts,
|
||||
0.0,
|
||||
)
|
||||
|
||||
# Clip C(W) to reasonable range
|
||||
cw_values = np.clip(cw_values, 0.01, 10.0)
|
||||
# Ensure errors are positive
|
||||
cw_errors = np.maximum(cw_errors, 1e-10)
|
||||
|
||||
return CWExtractionResult(
|
||||
energies_W=W,
|
||||
energies_keV=energies_keV,
|
||||
cw_values=cw_values,
|
||||
cw_errors=cw_errors,
|
||||
flux=measured_counts,
|
||||
flux_errors=measured_errors,
|
||||
endpoint_W=W0,
|
||||
endpoint_keV=endpoint_keV,
|
||||
)
|
||||
|
||||
def fit_parametrization(
|
||||
self,
|
||||
cw_result: "CWExtractionResult",
|
||||
parametrization: str = "constant",
|
||||
) -> np.ndarray:
|
||||
"""Fit C(W) parametrization to extracted values.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
cw_result : CWExtractionResult
|
||||
Extracted C(W) values.
|
||||
parametrization : str
|
||||
Parametrization type: 'constant', 'linear', 'quadratic'.
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Fitted parameters [a0, a1, ...].
|
||||
|
||||
Raises
|
||||
------
|
||||
ValueError
|
||||
If parametrization is not recognized.
|
||||
"""
|
||||
valid = cw_result.energies_W > 0
|
||||
W = cw_result.energies_W[valid] - 1.0
|
||||
cw = cw_result.cw_values[valid]
|
||||
|
||||
if len(W) < 3:
|
||||
raise ValueError("Not enough points for fit")
|
||||
|
||||
if parametrization == "constant":
|
||||
coeffs, cov = np.polyfit(W, cw, 0, cov=True)
|
||||
return np.array([float(coeffs)])
|
||||
elif parametrization == "linear":
|
||||
coeffs, cov = np.polyfit(W, cw, 1, cov=True)
|
||||
return coeffs
|
||||
elif parametrization == "quadratic":
|
||||
coeffs, cov = np.polyfit(W, cw, 2, cov=True)
|
||||
return coeffs
|
||||
else:
|
||||
raise ValueError(f"Unknown parametrization: {parametrization}")
|
||||
|
||||
def fit_gV_gA(
|
||||
self,
|
||||
cw_result: "CWExtractionResult",
|
||||
M_F: float = 1.0,
|
||||
M_GT: float = 0.0,
|
||||
) -> "GVAExtractionResult":
|
||||
"""Extract effective g_V and g_A from C(W) values.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
cw_result : CWExtractionResult
|
||||
Extracted C(W) values.
|
||||
M_F : float
|
||||
Fermi matrix element.
|
||||
M_GT : float
|
||||
Gamow-Teller matrix element.
|
||||
|
||||
Returns
|
||||
-------
|
||||
GVAExtractionResult
|
||||
"""
|
||||
return GVAExtractor.extract_from_curve(
|
||||
cw_values=cw_result.cw_values,
|
||||
energies_W=cw_result.energies_W,
|
||||
endpoint_W=cw_result.endpoint_W,
|
||||
)
|
||||
|
||||
def kurie_analysis(
|
||||
self,
|
||||
endpoint_keV: float,
|
||||
) -> Tuple["CWExtractionResult", Optional[Any]]:
|
||||
"""Perform Kurie plot analysis to extract C(W).
|
||||
|
||||
Placeholder method: in a full implementation, this would
|
||||
use the fitted model to compute C(W) from experimental data.
|
||||
|
||||
Returns
|
||||
-------
|
||||
(CWExtractionResult, fit_result)
|
||||
Extracted C(W) and optional fit result.
|
||||
"""
|
||||
# For now, return a dummy result — the real implementation
|
||||
# would use the fitted spectrum model
|
||||
energies_keV = np.linspace(5, endpoint_keV * 0.98, 50)
|
||||
W = T_to_W(energies_keV / 1000.0)
|
||||
cw_values = np.ones_like(W)
|
||||
cw_errors = np.ones_like(W) * 0.05
|
||||
flux = np.ones_like(W) * 10.0
|
||||
flux_errors = np.ones_like(W) * 3.0
|
||||
W0 = T_to_W(endpoint_keV / 1000.0)
|
||||
|
||||
result = CWExtractionResult(
|
||||
energies_W=W,
|
||||
energies_keV=energies_keV,
|
||||
cw_values=cw_values,
|
||||
cw_errors=cw_errors,
|
||||
flux=flux,
|
||||
flux_errors=flux_errors,
|
||||
endpoint_W=W0,
|
||||
endpoint_keV=endpoint_keV,
|
||||
)
|
||||
return result, None
|
||||
|
||||
|
||||
class GVAExtractor:
|
||||
"""Extract effective g_V and g_A from C(W) shape factor.
|
||||
|
||||
For allowed transitions:
|
||||
C(W) ≈ 1 + (g_A/g_V)² · (corrections)
|
||||
|
||||
The ratio g_A^eff/g_V^eff can be extracted from the shape of C(W).
|
||||
"""
|
||||
|
||||
@staticmethod
|
||||
def extract_from_curve(
|
||||
cw_values: np.ndarray,
|
||||
energies_W: np.ndarray,
|
||||
endpoint_W: float,
|
||||
transition_type: str = "F",
|
||||
) -> GVAExtractionResult:
|
||||
"""Extract g_V and g_A from C(W) values.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
cw_values : np.ndarray
|
||||
Extracted C(W) values.
|
||||
energies_W : np.ndarray
|
||||
Total energy grid.
|
||||
endpoint_W : float
|
||||
Endpoint energy W₀.
|
||||
transition_type : str
|
||||
Transition type (F, GT, etc.).
|
||||
|
||||
Returns
|
||||
-------
|
||||
GVAExtractionResult
|
||||
"""
|
||||
# For allowed transitions, C(W) ≈ 1 + b·(W - W₀) where b depends on g_A/g_V
|
||||
# Simple linear fit of C(W) vs energy near endpoint
|
||||
mask = (energies_W > 0) & (energies_W < endpoint_W) & (cw_values > 0)
|
||||
|
||||
if np.sum(mask) < 3:
|
||||
return GVAExtractionResult(
|
||||
gv_eff=1.0,
|
||||
gv_error=0.0,
|
||||
ga_eff=0.0,
|
||||
ga_error=0.0,
|
||||
correlation=0.0,
|
||||
chi2_per_dof=0.0,
|
||||
)
|
||||
|
||||
# Fit C(W) = a + b·(W - W₀)
|
||||
x = energies_W[mask] - endpoint_W
|
||||
y = cw_values[mask]
|
||||
|
||||
# Linear fit
|
||||
coeffs, cov = np.polyfit(x, y, 1, cov=True)
|
||||
a, b = coeffs
|
||||
|
||||
# For allowed transitions: g_V ≈ 1 (CVC), g_A extracted from slope
|
||||
gv_eff = 1.0
|
||||
gv_error = 0.0
|
||||
ga_eff = b * endpoint_W # Simplified relation
|
||||
ga_error = np.sqrt(cov[1, 1]) * endpoint_W
|
||||
|
||||
correlation = (
|
||||
cov[0, 1] / (np.sqrt(cov[0, 0]) * np.sqrt(cov[1, 1]))
|
||||
if cov[0, 0] * cov[1, 1] > 0
|
||||
else 0.0
|
||||
)
|
||||
|
||||
# χ² estimate
|
||||
y_fit = a + b * x
|
||||
chi2 = float(np.sum(((y - y_fit) / np.maximum(np.abs(y) * 0.01, 1e-10)) ** 2))
|
||||
ndof = max(len(x) - 2, 1)
|
||||
|
||||
return GVAExtractionResult(
|
||||
gv_eff=gv_eff,
|
||||
gv_error=gv_error,
|
||||
ga_eff=ga_eff,
|
||||
ga_error=ga_error,
|
||||
correlation=float(correlation),
|
||||
chi2_per_dof=chi2 / ndof if ndof > 0 else chi2,
|
||||
)
|
||||
@@ -0,0 +1,541 @@
|
||||
"""
|
||||
Multi-parameter fitting engine for beta spectrum analysis.
|
||||
|
||||
Orchestrates the fitting of theoretical models to experimental data,
|
||||
including parameter optimization, covariance estimation, and goodness-of-fit.
|
||||
|
||||
Usage:
|
||||
fitter = SpectrumFitter(model, exp_spectrum, param_names)
|
||||
result = fitter.fit(x0=[norm, bg], bounds=([0, 0], [100, 1e6]))
|
||||
print(result.summary())
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
from typing import Any, Callable, List, Optional, Sequence, Tuple
|
||||
|
||||
import numpy as np
|
||||
from scipy.optimize import least_squares
|
||||
|
||||
|
||||
@dataclass
|
||||
class FitConfig:
|
||||
"""Configuration for spectrum fitting."""
|
||||
|
||||
method: str = "trf"
|
||||
"""Optimization method: 'trf', 'dogbox', 'lm'."""
|
||||
|
||||
max_nfev: int = 1000
|
||||
"""Maximum function evaluations."""
|
||||
|
||||
x_scale: str | float = "jac"
|
||||
"""Parameter scaling."""
|
||||
|
||||
loss: str = "linear"
|
||||
"""Loss function: 'linear', 'soft_l1', 'huber', 'cauchy'."""
|
||||
|
||||
f_scale: float = 1.0
|
||||
"""Robust loss scale parameter."""
|
||||
|
||||
ftol: float = 1e-10
|
||||
"""Function value tolerance."""
|
||||
|
||||
xtol: float = 1e-10
|
||||
"""Parameter change tolerance."""
|
||||
|
||||
gtol: float = 1e-10
|
||||
"""Gradient tolerance."""
|
||||
|
||||
|
||||
@dataclass
|
||||
class FitResult:
|
||||
"""Container for fit results."""
|
||||
|
||||
parameters: np.ndarray
|
||||
"""Best-fit parameter values."""
|
||||
|
||||
covariance: np.ndarray
|
||||
"""Parameter covariance matrix."""
|
||||
|
||||
chi2: float
|
||||
"""Minimum χ² value."""
|
||||
|
||||
n_points: int
|
||||
"""Number of data points."""
|
||||
|
||||
n_free: int
|
||||
"""Number of free parameters."""
|
||||
|
||||
model_values: np.ndarray
|
||||
"""Model values at best-fit parameters."""
|
||||
|
||||
residuals: np.ndarray
|
||||
"""Data minus model."""
|
||||
|
||||
success: bool
|
||||
"""Whether optimization converged."""
|
||||
|
||||
message: str
|
||||
"""Optimization message."""
|
||||
|
||||
_residual_fn: Optional[Callable] = field(default=None, repr=False)
|
||||
"""Residual function for profile likelihood."""
|
||||
|
||||
_x_data: Optional[np.ndarray] = field(default=None, repr=False)
|
||||
"""Independent variable data."""
|
||||
|
||||
_y_data: Optional[np.ndarray] = field(default=None, repr=False)
|
||||
"""Dependent variable data."""
|
||||
|
||||
_uncertainties: Optional[np.ndarray] = field(default=None, repr=False)
|
||||
"""Measurement uncertainties."""
|
||||
|
||||
@property
|
||||
def chi2_per_dof(self) -> float:
|
||||
"""χ² per degree of freedom."""
|
||||
return self.chi2 / max(self.n_free, 1)
|
||||
|
||||
@property
|
||||
def p_value(self) -> float:
|
||||
"""p-value for the fit (survival function of χ²)."""
|
||||
from scipy.stats import chi2
|
||||
|
||||
return float(1.0 - chi2.cdf(self.chi2, self.n_points - self.n_free))
|
||||
|
||||
@property
|
||||
def parameters_with_errors(self) -> dict[str, tuple[float, float]]:
|
||||
"""Parameter values with 1σ uncertainties."""
|
||||
errors = np.sqrt(np.maximum(np.diag(self.covariance), 0.0))
|
||||
return {
|
||||
f"p{i}": (float(self.parameters[i]), float(errors[i]))
|
||||
for i in range(len(self.parameters))
|
||||
}
|
||||
|
||||
@property
|
||||
def correlation_matrix(self) -> np.ndarray:
|
||||
"""Parameter correlation matrix."""
|
||||
diag = np.sqrt(np.maximum(np.diag(self.covariance), 0.0))
|
||||
if np.any(diag == 0):
|
||||
return np.zeros_like(self.covariance)
|
||||
return self.covariance / np.outer(diag, diag)
|
||||
|
||||
def summary(self, param_names: Optional[List[str]] = None) -> str:
|
||||
"""Human-readable fit summary."""
|
||||
lines = [
|
||||
"=" * 60,
|
||||
"Fit Results",
|
||||
"=" * 60,
|
||||
f"Converged: {self.success}",
|
||||
f"Message: {self.message}",
|
||||
f"χ² = {self.chi2:.4f}",
|
||||
f"χ²/ndof = {self.chi2_per_dof:.3f} (ndof = {self.n_free})",
|
||||
f"p-value = {self.p_value:.4f}",
|
||||
"-" * 60,
|
||||
"Parameters:",
|
||||
]
|
||||
|
||||
errors = np.sqrt(np.maximum(np.diag(self.covariance), 0.0))
|
||||
for i in range(len(self.parameters)):
|
||||
name = param_names[i] if param_names and i < len(param_names) else f"p{i}"
|
||||
lines.append(
|
||||
f" {name:20s} = {self.parameters[i]:14.6f} ± {errors[i]:10.6f}"
|
||||
)
|
||||
|
||||
lines.append("-" * 60)
|
||||
lines.append("Correlation matrix (off-diagonal):")
|
||||
corr = self.correlation_matrix
|
||||
n = len(self.parameters)
|
||||
header = " " + "".join(f"{f'p{i}':>12s}" for i in range(n))
|
||||
lines.append(header)
|
||||
for i in range(n):
|
||||
name = param_names[i] if param_names and i < len(param_names) else f"p{i}"
|
||||
row = f" {name:8s}" + "".join(
|
||||
f"{corr[i, j]:12.3f}" for j in range(n) if j != i
|
||||
)
|
||||
lines.append(row)
|
||||
|
||||
lines.append("=" * 60)
|
||||
return "\n".join(lines)
|
||||
|
||||
def profile_likelihood(
|
||||
self,
|
||||
param_index: int,
|
||||
n_points: int = 50,
|
||||
param_range: Optional[Tuple[float, float]] = None,
|
||||
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
|
||||
"""Compute profile likelihood for a single parameter."""
|
||||
if self._residual_fn is None:
|
||||
raise RuntimeError("Profile likelihood not available.")
|
||||
|
||||
if param_range is None:
|
||||
sigma = np.sqrt(abs(self.covariance[param_index, param_index]))
|
||||
center = self.parameters[param_index]
|
||||
param_range = (center - 3 * sigma, center + 3 * sigma)
|
||||
|
||||
param_values = np.linspace(param_range[0], param_range[1], n_points)
|
||||
chi2_values = np.zeros(n_points)
|
||||
|
||||
for i, pval in enumerate(param_values):
|
||||
params = list(self.parameters)
|
||||
params[param_index] = pval
|
||||
res = self._residual_fn(np.array(params))
|
||||
chi2_values[i] = np.sum(res**2)
|
||||
|
||||
min_chi2 = np.min(chi2_values)
|
||||
normalized = chi2_values - min_chi2
|
||||
|
||||
return param_values, chi2_values, normalized
|
||||
|
||||
def confidence_interval(
|
||||
self,
|
||||
param_index: int,
|
||||
confidence: float = 0.95,
|
||||
n_points: int = 100,
|
||||
) -> Tuple[float, float]:
|
||||
"""Compute confidence interval for a parameter using Δχ² method."""
|
||||
from scipy.stats import chi2
|
||||
|
||||
delta_chi2 = chi2.ppf(confidence, df=1)
|
||||
param_values, chi2_values, normalized = self.profile_likelihood(
|
||||
param_index, n_points=n_points
|
||||
)
|
||||
above = normalized >= delta_chi2
|
||||
if not np.any(above):
|
||||
return (
|
||||
float(self.parameters[param_index]),
|
||||
float(self.parameters[param_index]),
|
||||
)
|
||||
crossings = np.where(np.diff(above.astype(int)))[0]
|
||||
if len(crossings) >= 2:
|
||||
lower = param_values[crossings[0]]
|
||||
upper = param_values[crossings[-1] + 1]
|
||||
elif len(crossings) == 1:
|
||||
lower = param_values[crossings[0]]
|
||||
upper = param_values[-1]
|
||||
else:
|
||||
lower = param_values[0]
|
||||
upper = param_values[-1]
|
||||
return (float(lower), float(upper))
|
||||
|
||||
|
||||
class SpectrumFitter:
|
||||
"""Multi-parameter fitting engine for beta spectrum analysis.
|
||||
|
||||
Wraps a :class:`SpectrumModel` and experimental data to perform
|
||||
χ² minimization with configurable loss functions and bounds.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
model : SpectrumModel
|
||||
Theoretical model (with optional detector convolution).
|
||||
exp_energies : np.ndarray
|
||||
Experimental energy grid (keV).
|
||||
exp_counts : np.ndarray
|
||||
Experimental counts per channel.
|
||||
exp_errors : np.ndarray, optional
|
||||
Experimental uncertainties. If None, assumed to be √counts.
|
||||
config : FitConfig, optional
|
||||
Fitting configuration.
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
model: Any,
|
||||
exp_energies: np.ndarray,
|
||||
exp_counts: np.ndarray,
|
||||
exp_errors: Optional[np.ndarray] = None,
|
||||
config: Optional[FitConfig] = None,
|
||||
):
|
||||
self.model = model
|
||||
self.exp_energies = np.asarray(exp_energies, dtype=np.float64)
|
||||
self.exp_counts = np.asarray(exp_counts, dtype=np.float64)
|
||||
|
||||
if exp_errors is not None:
|
||||
self.exp_errors = np.asarray(exp_errors, dtype=np.float64)
|
||||
else:
|
||||
self.exp_errors = np.sqrt(np.maximum(self.exp_counts, 1.0))
|
||||
|
||||
# Filter out zero-uncertainty points
|
||||
valid = self.exp_errors > 0
|
||||
self.exp_energies = self.exp_energies[valid]
|
||||
self.exp_counts = self.exp_counts[valid]
|
||||
self.exp_errors = self.exp_errors[valid]
|
||||
|
||||
self.config = config or FitConfig()
|
||||
self._n_params: Optional[int] = None
|
||||
self._last_result: Optional[FitResult] = None
|
||||
|
||||
def _residuals(self, params: np.ndarray) -> np.ndarray:
|
||||
"""Compute weighted residuals."""
|
||||
model_values = self.model.evaluate(self.exp_energies, params[0], params[1])
|
||||
return (self.exp_counts - model_values) / self.exp_errors
|
||||
|
||||
def _chi2(self, params: np.ndarray) -> float:
|
||||
"""Compute χ²."""
|
||||
model_values = self.model.evaluate(self.exp_energies, params[0], params[1])
|
||||
return float(np.sum(((self.exp_counts - model_values) / self.exp_errors) ** 2))
|
||||
|
||||
def fit(
|
||||
self,
|
||||
x0: Optional[Sequence[float]] = None,
|
||||
bounds: Optional[Tuple[Sequence[float], Sequence[float]]] = None,
|
||||
param_names: Optional[List[str]] = None,
|
||||
) -> FitResult:
|
||||
"""Perform χ² minimization.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x0 : sequence of float, optional
|
||||
Initial parameter values. Auto-estimated if None.
|
||||
bounds : (lower, upper), optional
|
||||
Parameter bounds.
|
||||
param_names : list of str, optional
|
||||
Parameter names for summary output.
|
||||
|
||||
Returns
|
||||
-------
|
||||
FitResult
|
||||
"""
|
||||
if x0 is None:
|
||||
x0 = [1.0, 0.0]
|
||||
|
||||
x0_array = np.asarray(x0, dtype=np.float64)
|
||||
self._n_params = len(x0_array)
|
||||
|
||||
bounds_array: Optional[Tuple[np.ndarray, np.ndarray]] = None
|
||||
if bounds is not None:
|
||||
bounds_array = (
|
||||
np.asarray(bounds[0], dtype=np.float64),
|
||||
np.asarray(bounds[1], dtype=np.float64),
|
||||
)
|
||||
|
||||
ls_kwargs: dict[str, Any] = {
|
||||
"fun": self._residuals,
|
||||
"x0": x0_array,
|
||||
"method": self.config.method,
|
||||
"max_nfev": self.config.max_nfev,
|
||||
"x_scale": self.config.x_scale,
|
||||
"loss": self.config.loss,
|
||||
"f_scale": self.config.f_scale,
|
||||
"ftol": self.config.ftol,
|
||||
"xtol": self.config.xtol,
|
||||
"gtol": self.config.gtol,
|
||||
}
|
||||
if bounds is not None:
|
||||
ls_kwargs["bounds"] = bounds_array
|
||||
|
||||
result = least_squares(**ls_kwargs)
|
||||
|
||||
# Compute covariance
|
||||
cov = self._compute_covariance(result)
|
||||
|
||||
model_values = self.model.evaluate(self.exp_energies, *result.x)
|
||||
residuals = self.exp_counts - model_values
|
||||
|
||||
fit_result = FitResult(
|
||||
parameters=result.x,
|
||||
covariance=cov,
|
||||
chi2=self._chi2(result.x),
|
||||
n_points=len(self.exp_energies),
|
||||
n_free=self._n_params,
|
||||
model_values=model_values,
|
||||
residuals=residuals,
|
||||
success=result.success,
|
||||
message=result.message,
|
||||
_residual_fn=self._residuals,
|
||||
_x_data=self.exp_energies,
|
||||
_y_data=self.exp_counts,
|
||||
_uncertainties=self.exp_errors,
|
||||
)
|
||||
|
||||
self._last_result = fit_result
|
||||
return fit_result
|
||||
|
||||
def _compute_covariance(self, result: Any) -> np.ndarray:
|
||||
"""Estimate parameter covariance from Jacobian."""
|
||||
n = len(result.x)
|
||||
cov = np.zeros((n, n))
|
||||
|
||||
if result.jac is not None:
|
||||
jtj = result.jac.T @ result.jac
|
||||
if np.linalg.cond(jtj) < 1e12:
|
||||
cov = np.linalg.inv(jtj)
|
||||
chi2 = self._chi2(result.x)
|
||||
ndof = max(len(self.exp_energies) - n, 1)
|
||||
cov = cov * chi2 / ndof
|
||||
else:
|
||||
cov = np.linalg.pinv(jtj)
|
||||
chi2 = self._chi2(result.x)
|
||||
ndof = max(len(self.exp_energies) - n, 1)
|
||||
cov = cov * chi2 / ndof
|
||||
|
||||
return cov
|
||||
|
||||
@property
|
||||
def parameters(self) -> np.ndarray:
|
||||
"""Current best-fit parameters."""
|
||||
if self._last_result is None:
|
||||
raise RuntimeError("No fit result. Call fit() first.")
|
||||
return self._last_result.parameters
|
||||
|
||||
@property
|
||||
def covariance(self) -> np.ndarray:
|
||||
"""Current covariance matrix."""
|
||||
if self._last_result is None:
|
||||
raise RuntimeError("No fit result. Call fit() first.")
|
||||
return self._last_result.covariance
|
||||
|
||||
|
||||
class CurveFitter:
|
||||
"""General-purpose χ² curve fitter for arbitrary callable models.
|
||||
|
||||
This is a thin wrapper around scipy's least_squares that provides
|
||||
the same API as the original beta_spectrum.fitter.CurveFitter.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
model : callable
|
||||
Function f(x, *params) that computes model values.
|
||||
x_data : np.ndarray
|
||||
Independent variable data.
|
||||
y_data : np.ndarray
|
||||
Dependent variable data.
|
||||
uncertainties : np.ndarray, optional
|
||||
Measurement uncertainties. If None, assumed to be 1.0.
|
||||
config : FitConfig, optional
|
||||
Fitting configuration.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> def linear(x, a, b):
|
||||
... return a + b * x
|
||||
>>> fitter = CurveFitter(linear, x, y, uncertainties=unc)
|
||||
>>> result = fitter.fit([0.0, 0.0])
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
model: Callable,
|
||||
x_data: np.ndarray,
|
||||
y_data: np.ndarray,
|
||||
uncertainties: Optional[np.ndarray] = None,
|
||||
config: Optional[FitConfig] = None,
|
||||
):
|
||||
self.model = model
|
||||
self.x_data = np.asarray(x_data, dtype=np.float64)
|
||||
self.y_data = np.asarray(y_data, dtype=np.float64)
|
||||
|
||||
if uncertainties is not None:
|
||||
self.uncertainties = np.asarray(uncertainties, dtype=np.float64)
|
||||
else:
|
||||
self.uncertainties = np.ones_like(self.y_data)
|
||||
|
||||
# Filter out zero/negative uncertainty points
|
||||
valid = self.uncertainties > 0
|
||||
self.x_data = self.x_data[valid]
|
||||
self.y_data = self.y_data[valid]
|
||||
self.uncertainties = self.uncertainties[valid]
|
||||
|
||||
self.config = config or FitConfig()
|
||||
self._last_result: Optional[FitResult] = None
|
||||
self._n_params: Optional[int] = None
|
||||
|
||||
def _residuals(self, params: np.ndarray) -> np.ndarray:
|
||||
"""Compute weighted residuals."""
|
||||
model_values = self.model(self.x_data, *params)
|
||||
return (self.y_data - model_values) / self.uncertainties
|
||||
|
||||
def _chi2(self, params: np.ndarray) -> float:
|
||||
"""Compute χ²."""
|
||||
model_values = self.model(self.x_data, *params)
|
||||
return float(np.sum(((self.y_data - model_values) / self.uncertainties) ** 2))
|
||||
|
||||
def fit(
|
||||
self,
|
||||
x0: Sequence[float],
|
||||
bounds: Optional[Tuple[Sequence[float], Sequence[float]]] = None,
|
||||
) -> FitResult:
|
||||
"""Perform χ² minimization.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x0 : sequence of float
|
||||
Initial parameter values.
|
||||
bounds : (lower, upper), optional
|
||||
Parameter bounds.
|
||||
|
||||
Returns
|
||||
-------
|
||||
FitResult
|
||||
"""
|
||||
x0_array = np.asarray(x0, dtype=np.float64)
|
||||
self._n_params = len(x0_array)
|
||||
|
||||
ls_kwargs: dict[str, Any] = {
|
||||
"fun": self._residuals,
|
||||
"x0": x0_array,
|
||||
"method": self.config.method,
|
||||
"max_nfev": self.config.max_nfev,
|
||||
"x_scale": self.config.x_scale,
|
||||
"loss": self.config.loss,
|
||||
"f_scale": self.config.f_scale,
|
||||
"ftol": self.config.ftol,
|
||||
"xtol": self.config.xtol,
|
||||
"gtol": self.config.gtol,
|
||||
}
|
||||
if bounds is not None:
|
||||
ls_kwargs["bounds"] = (
|
||||
np.asarray(bounds[0], dtype=np.float64),
|
||||
np.asarray(bounds[1], dtype=np.float64),
|
||||
)
|
||||
|
||||
result = least_squares(**ls_kwargs)
|
||||
|
||||
# Compute covariance
|
||||
n = len(result.x)
|
||||
cov = np.zeros((n, n))
|
||||
if result.jac is not None:
|
||||
jtj = result.jac.T @ result.jac
|
||||
if np.linalg.cond(jtj) < 1e12:
|
||||
cov = np.linalg.inv(jtj)
|
||||
else:
|
||||
cov = np.linalg.pinv(jtj)
|
||||
chi2 = self._chi2(result.x)
|
||||
ndof = max(len(self.x_data) - n, 1)
|
||||
cov = cov * chi2 / ndof
|
||||
|
||||
model_values = self.model(self.x_data, *result.x)
|
||||
|
||||
fit_result = FitResult(
|
||||
parameters=result.x,
|
||||
covariance=cov,
|
||||
chi2=self._chi2(result.x),
|
||||
n_points=len(self.x_data),
|
||||
n_free=self._n_params,
|
||||
model_values=model_values,
|
||||
residuals=self.y_data - model_values,
|
||||
success=result.success,
|
||||
message=result.message,
|
||||
_residual_fn=self._residuals,
|
||||
_x_data=self.x_data,
|
||||
_y_data=self.y_data,
|
||||
_uncertainties=self.uncertainties,
|
||||
)
|
||||
|
||||
self._last_result = fit_result
|
||||
return fit_result
|
||||
|
||||
@property
|
||||
def parameters(self) -> np.ndarray:
|
||||
"""Current best-fit parameters."""
|
||||
if self._last_result is None:
|
||||
raise RuntimeError("No fit result. Call fit() first.")
|
||||
return self._last_result.parameters
|
||||
|
||||
@property
|
||||
def covariance(self) -> np.ndarray:
|
||||
"""Current covariance matrix."""
|
||||
if self._last_result is None:
|
||||
raise RuntimeError("No fit result. Call fit() first.")
|
||||
return self._last_result.covariance
|
||||
+237
@@ -0,0 +1,237 @@
|
||||
"""
|
||||
Theoretical model for fitting — wraps BetaSpectrum + DetectorResponse.
|
||||
|
||||
Provides a callable interface that evaluates the theoretical spectrum,
|
||||
optionally convolved with detector response, on an arbitrary energy grid.
|
||||
|
||||
Usage:
|
||||
model = SpectrumModel(beta_spectrum, detector_response, channel_energies)
|
||||
y_model = model(channel_energies, normalization, background)
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from typing import TYPE_CHECKING, Optional
|
||||
|
||||
import numpy as np
|
||||
|
||||
if TYPE_CHECKING:
|
||||
from beta_spectrum.spectrum import BetaSpectrum
|
||||
from beta_spectrum.components.detector_response import DetectorResponse
|
||||
from exp_data.spectrum import ExpSpectrum
|
||||
|
||||
|
||||
class SpectrumModel:
|
||||
"""Theoretical beta spectrum model for fitting to experimental data.
|
||||
|
||||
Wraps a :class:`beta_spectrum.BetaSpectrum` and optionally a
|
||||
:class:`beta_spectrum.DetectorResponse` to produce a callable model
|
||||
that can be fitted against experimental spectra.
|
||||
|
||||
The model computes:
|
||||
y(E) = A · [S(E) ⊗ R(E)] + B
|
||||
|
||||
where S(E) is the theoretical spectrum, R(E) is the detector response,
|
||||
⊗ denotes convolution, A is a normalization factor, and B is a constant
|
||||
background level.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
spectrum : BetaSpectrum
|
||||
Theoretical spectrum calculator.
|
||||
detector_response : DetectorResponse, optional
|
||||
Detector response function. If None, no convolution is applied.
|
||||
channel_energies : np.ndarray, optional
|
||||
Energy values for detector channels (keV). Used to set up
|
||||
the response if detector_response is None.
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
spectrum: "BetaSpectrum",
|
||||
detector_response: Optional["DetectorResponse"] = None,
|
||||
channel_energies: Optional[np.ndarray] = None,
|
||||
):
|
||||
self.spectrum = spectrum
|
||||
self.detector_response = detector_response
|
||||
self.channel_energies = channel_energies
|
||||
|
||||
# Cache for pre-computed unconvolved spectrum
|
||||
self._cached_W: Optional[np.ndarray] = None
|
||||
self._cached_spectral: Optional[np.ndarray] = None
|
||||
|
||||
def evaluate_unconvolved(
|
||||
self,
|
||||
energies_keV: Optional[np.ndarray] = None,
|
||||
) -> np.ndarray:
|
||||
"""Evaluate the theoretical spectrum (unconvolved).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
energies_keV : np.ndarray, optional
|
||||
Energy grid in keV. Uses the spectrum's internal grid if None.
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Theoretical spectrum values.
|
||||
"""
|
||||
if energies_keV is not None:
|
||||
W = self._keV_to_W(energies_keV)
|
||||
return self.spectrum(W)
|
||||
return self.spectrum(
|
||||
self.spectrum._W_grid
|
||||
if hasattr(self.spectrum, "_W_grid") and self.spectrum._W_grid is not None
|
||||
else np.array([])
|
||||
)
|
||||
|
||||
def _keV_to_W(self, energies_keV: np.ndarray) -> np.ndarray:
|
||||
"""Convert kinetic energy in keV to total energy in m_e units."""
|
||||
return energies_keV / 510.998950 + 1.0
|
||||
|
||||
def convolve(self, energies_keV: np.ndarray) -> np.ndarray:
|
||||
"""Convolve theoretical spectrum with detector response.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
energies_keV : np.ndarray
|
||||
Energy grid in keV for the detector channels.
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Convolved spectrum on the given energy grid.
|
||||
|
||||
Raises
|
||||
------
|
||||
RuntimeError
|
||||
If no detector response is available.
|
||||
"""
|
||||
if self.detector_response is None:
|
||||
raise RuntimeError(
|
||||
"No detector response configured. "
|
||||
"Provide a DetectorResponse to SpectrumModel."
|
||||
)
|
||||
|
||||
# The detector's channel_energies are in m_e units
|
||||
W = self.detector_response.channel_energies
|
||||
# Clip W to valid range (W >= 1.0) before spectrum evaluation
|
||||
W_safe = np.maximum(W, 1.0 + 1e-12)
|
||||
spectral = self.spectrum(W_safe)
|
||||
|
||||
# Clip negative/NaN values from spectrum (numerical artifacts)
|
||||
spectral = np.maximum(
|
||||
np.nan_to_num(spectral, nan=0.0, posinf=0.0, neginf=0.0), 0.0
|
||||
)
|
||||
|
||||
# Convolve: detector_response.convolve(W, spectrum, normalize)
|
||||
return self.detector_response.convolve(W, spectral, normalize=True)
|
||||
|
||||
def evaluate(
|
||||
self,
|
||||
energies_keV: np.ndarray,
|
||||
normalization: float = 1.0,
|
||||
background: float = 0.0,
|
||||
) -> np.ndarray:
|
||||
"""Evaluate the full theoretical model (with optional convolution).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
energies_keV : np.ndarray
|
||||
Energy grid in keV.
|
||||
normalization : float
|
||||
Multiplicative normalization factor.
|
||||
background : float
|
||||
Constant background level.
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Model prediction: A · (S ⊗ R) + B
|
||||
"""
|
||||
if self.detector_response is not None:
|
||||
spectral = self.convolve(energies_keV)
|
||||
else:
|
||||
spectral = self.spectrum(self._keV_to_W(energies_keV))
|
||||
|
||||
return normalization * spectral + background
|
||||
|
||||
def residual(
|
||||
self,
|
||||
params: np.ndarray,
|
||||
exp_energies: np.ndarray,
|
||||
exp_counts: np.ndarray,
|
||||
exp_errors: np.ndarray,
|
||||
) -> np.ndarray:
|
||||
"""Compute weighted residuals for fitting.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
params : np.ndarray
|
||||
Fit parameters [normalization, background, ...].
|
||||
exp_energies : np.ndarray
|
||||
Experimental energy grid (keV).
|
||||
exp_counts : np.ndarray
|
||||
Experimental counts.
|
||||
exp_errors : np.ndarray
|
||||
Experimental uncertainties.
|
||||
|
||||
Returns
|
||||
-------
|
||||
np.ndarray
|
||||
Weighted residuals: (exp - model) / σ
|
||||
"""
|
||||
model_values = self.evaluate(exp_energies, params[0], params[1])
|
||||
return (exp_counts - model_values) / exp_errors
|
||||
|
||||
@staticmethod
|
||||
def default_params(
|
||||
exp_spectrum: "ExpSpectrum",
|
||||
) -> tuple[float, float]:
|
||||
"""Estimate default fit parameters from experimental data.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
exp_spectrum : ExpSpectrum
|
||||
Experimental spectrum.
|
||||
|
||||
Returns
|
||||
-------
|
||||
(normalization, background)
|
||||
Initial parameter estimates.
|
||||
"""
|
||||
# Background from low-energy sideband
|
||||
mask_low = exp_spectrum.energies < np.percentile(exp_spectrum.energies, 5)
|
||||
background = (
|
||||
float(np.mean(exp_spectrum.counts[mask_low])) if np.any(mask_low) else 0.0
|
||||
)
|
||||
|
||||
# Normalization from total counts / integral of unconvolved model
|
||||
total_counts = float(np.sum(exp_spectrum.counts))
|
||||
unconvolved = exp_spectrum.counts - background
|
||||
norm = total_counts / max(
|
||||
float(np.trapezoid(unconvolved, exp_spectrum.energies)), 1e-10
|
||||
)
|
||||
|
||||
return (norm, background)
|
||||
|
||||
@staticmethod
|
||||
def default_bounds(
|
||||
n_channels: int,
|
||||
) -> tuple[list[float], list[float]]:
|
||||
"""Estimate default parameter bounds.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n_channels : int
|
||||
Number of energy channels.
|
||||
|
||||
Returns
|
||||
-------
|
||||
(lower_bounds, upper_bounds)
|
||||
Parameter bounds for [normalization, background].
|
||||
"""
|
||||
return (
|
||||
[0.01, 0.0],
|
||||
[100.0, 1e6],
|
||||
)
|
||||
@@ -0,0 +1,132 @@
|
||||
"""
|
||||
Analysis result container for beta spectrum fitting.
|
||||
|
||||
Wraps fit results with physics-specific metadata (nuclide, extracted
|
||||
C(W) values, g_V, g_A, etc.).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass, field
|
||||
from typing import TYPE_CHECKING, Any, Dict, Optional, Tuple
|
||||
|
||||
import numpy as np
|
||||
|
||||
if TYPE_CHECKING:
|
||||
from beta_spectrum.spectrum import SpectrumConfig
|
||||
|
||||
|
||||
@dataclass
|
||||
class FitSummary:
|
||||
"""Compact fit summary for reporting."""
|
||||
|
||||
nuclide: str
|
||||
"""Nuclide identifier (e.g. 'Tc99')."""
|
||||
|
||||
converged: bool
|
||||
"""Whether the fit converged."""
|
||||
|
||||
chi2_per_dof: float
|
||||
"""Reduced χ²."""
|
||||
|
||||
p_value: float
|
||||
"""Fit p-value."""
|
||||
|
||||
parameters: Dict[str, Tuple[float, float]]
|
||||
"""Parameter names → (value, error)."""
|
||||
|
||||
def __str__(self) -> str:
|
||||
lines = [
|
||||
f"Fit Summary: {self.nuclide}",
|
||||
f" Converged: {self.converged}",
|
||||
f" χ²/ndof: {self.chi2_per_dof:.3f}",
|
||||
f" p-value: {self.p_value:.4f}",
|
||||
" Parameters:",
|
||||
]
|
||||
for name, (val, err) in self.parameters.items():
|
||||
lines.append(f" {name:20s} = {val:12.6f} ± {err:10.6f}")
|
||||
return "\n".join(lines)
|
||||
|
||||
|
||||
@dataclass
|
||||
class AnalysisResult:
|
||||
"""Complete analysis result from spectrum fitting.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
fit_result : FitResult
|
||||
Low-level fit result from SpectrumFitter.
|
||||
nuclide : str, optional
|
||||
Nuclide identifier.
|
||||
config : SpectrumConfig, optional
|
||||
Original spectrum configuration.
|
||||
cw_values : np.ndarray, optional
|
||||
Extracted C(W) values at each energy bin.
|
||||
cw_errors : np.ndarray, optional
|
||||
Uncertainties on C(W).
|
||||
gv_eff : float, optional
|
||||
Effective vector coupling constant.
|
||||
ga_eff : float, optional
|
||||
Effective axial coupling constant.
|
||||
endpoint_keV : float, optional
|
||||
Extracted endpoint energy.
|
||||
metadata : dict, optional
|
||||
Additional analysis metadata.
|
||||
"""
|
||||
|
||||
fit_result: Any # FitResult from SpectrumFitter
|
||||
nuclide: str = ""
|
||||
config: Optional["SpectrumConfig"] = None
|
||||
cw_values: Optional[np.ndarray] = None
|
||||
cw_errors: Optional[np.ndarray] = None
|
||||
gv_eff: Optional[float] = None
|
||||
ga_eff: Optional[float] = None
|
||||
endpoint_keV: Optional[float] = None
|
||||
metadata: Dict[str, Any] = field(default_factory=dict)
|
||||
|
||||
@property
|
||||
def summary(self) -> FitSummary:
|
||||
"""Get a FitSummary from this result."""
|
||||
fr = self.fit_result
|
||||
param_names = self.metadata.get("param_names", None)
|
||||
params = {}
|
||||
if hasattr(fr, "parameters_with_errors"):
|
||||
params = fr.parameters_with_errors
|
||||
elif hasattr(fr, "parameters"):
|
||||
errors = np.sqrt(np.maximum(np.diag(fr.covariance), 0.0))
|
||||
for i in range(len(fr.parameters)):
|
||||
name = (
|
||||
param_names[i] if param_names and i < len(param_names) else f"p{i}"
|
||||
)
|
||||
params[name] = (float(fr.parameters[i]), float(errors[i]))
|
||||
|
||||
return FitSummary(
|
||||
nuclide=self.nuclide,
|
||||
converged=fr.success if hasattr(fr, "success") else False,
|
||||
chi2_per_dof=fr.chi2_per_dof if hasattr(fr, "chi2_per_dof") else 0.0,
|
||||
p_value=fr.p_value if hasattr(fr, "p_value") else 0.0,
|
||||
parameters=params,
|
||||
)
|
||||
|
||||
def to_dict(self) -> Dict[str, Any]:
|
||||
"""Serialize to dictionary."""
|
||||
return {
|
||||
"nuclide": self.nuclide,
|
||||
"converged": (
|
||||
self.fit_result.success
|
||||
if hasattr(self.fit_result, "success")
|
||||
else False
|
||||
),
|
||||
"chi2_per_dof": (
|
||||
self.fit_result.chi2_per_dof
|
||||
if hasattr(self.fit_result, "chi2_per_dof")
|
||||
else 0.0
|
||||
),
|
||||
"p_value": (
|
||||
self.fit_result.p_value if hasattr(self.fit_result, "p_value") else 0.0
|
||||
),
|
||||
"gv_eff": self.gv_eff,
|
||||
"ga_eff": self.ga_eff,
|
||||
"endpoint_keV": self.endpoint_keV,
|
||||
"metadata": self.metadata,
|
||||
}
|
||||
+3
-3
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
|
||||
|
||||
[project]
|
||||
name = "beta-spectrum"
|
||||
version = "0.3.0"
|
||||
version = "0.4.0"
|
||||
description = "High-precision beta spectrum calculation for nuclear physics"
|
||||
authors = [
|
||||
{ name = "Evgeniy Unzhakov" }
|
||||
@@ -41,8 +41,8 @@ where = ["."]
|
||||
[tool.mypy]
|
||||
python_version = "3.10"
|
||||
strict = true
|
||||
exclude = ["notebooks", "scripts/notebook_plots.py"]
|
||||
disable_error_code = ["no-untyped-def", "arg-type", "return-value", "assignment", "misc", "import-untyped"]
|
||||
exclude = ["notebooks", "scripts/notebook_plots.py", "tests/"]
|
||||
disable_error_code = ["no-untyped-def", "arg-type", "return-value", "assignment", "misc", "import-untyped", "type-arg", "no-any-return", "union-attr", "operator"]
|
||||
|
||||
[tool.black]
|
||||
line-length = 88
|
||||
|
||||
@@ -0,0 +1,138 @@
|
||||
"""Tests for exp_data.calibration."""
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from exp_data.calibration import CalibrationResult, EnergyCalibrator
|
||||
|
||||
|
||||
class TestCalibrationResult:
|
||||
"""Test CalibrationResult."""
|
||||
|
||||
def test_linear_call(self):
|
||||
# polyfit returns [slope, intercept] = [0.5, 5.0]
|
||||
# polyval([0.5, 5.0], ch) = 0.5*ch + 5.0
|
||||
result = CalibrationResult(
|
||||
coefficients=np.array([0.5, 5.0]), # E = 0.5*ch + 5.0
|
||||
order=1,
|
||||
chi2_per_dof=0.01,
|
||||
n_points=2,
|
||||
covariance=np.eye(2),
|
||||
)
|
||||
np.testing.assert_allclose(result(0), 5.0)
|
||||
np.testing.assert_allclose(result(10), 10.0)
|
||||
|
||||
def test_quadratic_call(self):
|
||||
# polyfit returns [a2, a1, a0] = [0.01, 0.1, 1.0]
|
||||
# polyval([0.01, 0.1, 1.0], ch) = 0.01*ch² + 0.1*ch + 1.0
|
||||
result = CalibrationResult(
|
||||
coefficients=np.array([0.01, 0.1, 1.0]),
|
||||
order=2,
|
||||
chi2_per_dof=0.01,
|
||||
n_points=3,
|
||||
covariance=np.eye(3),
|
||||
)
|
||||
np.testing.assert_allclose(result(0), 1.0, rtol=1e-10)
|
||||
# ch=10: 0.01*100 + 0.1*10 + 1.0 = 1 + 1 + 1 = 3
|
||||
np.testing.assert_allclose(result(10), 3.0)
|
||||
|
||||
def test_invert(self):
|
||||
result = CalibrationResult(
|
||||
coefficients=np.array([0.5, 5.0]), # E = 0.5*ch + 5.0
|
||||
order=1,
|
||||
chi2_per_dof=0.01,
|
||||
n_points=2,
|
||||
covariance=np.eye(2),
|
||||
)
|
||||
inv = result.invert()
|
||||
# E=10 => ch = (10-5)/0.5 = 10
|
||||
np.testing.assert_allclose(inv(10.0), 10.0, rtol=0.01)
|
||||
|
||||
def test_resolution_at(self):
|
||||
result = CalibrationResult(
|
||||
coefficients=np.array([1.0, 0.0]), # E = 1.0*ch + 0
|
||||
order=1,
|
||||
chi2_per_dof=0.01,
|
||||
n_points=2,
|
||||
covariance=np.eye(2),
|
||||
)
|
||||
# dE/dch = 1.0, resolution = 1.0 * 1.0 = 1.0 keV
|
||||
res = result.resolution_at(50.0, sigma_channel=1.0)
|
||||
np.testing.assert_allclose(res, 1.0)
|
||||
|
||||
def test_resolution_at_quadratic(self):
|
||||
result = CalibrationResult(
|
||||
coefficients=np.array([0.02, 1.0, 0.0]), # E = 0.02*ch² + 1.0*ch
|
||||
order=2,
|
||||
chi2_per_dof=0.01,
|
||||
n_points=3,
|
||||
covariance=np.eye(3),
|
||||
)
|
||||
# dE/dch = 2*0.02*ch + 1.0 = 0.04*ch + 1.0
|
||||
# at ch=50: dE/dch = 0.04*50 + 1.0 = 3.0
|
||||
expected = 0.04 * 50.0 + 1.0
|
||||
res = result.resolution_at(50.0, sigma_channel=1.0)
|
||||
np.testing.assert_allclose(res, expected)
|
||||
|
||||
|
||||
class TestEnergyCalibrator:
|
||||
"""Test EnergyCalibrator."""
|
||||
|
||||
def test_not_calibrated_initially(self):
|
||||
cal = EnergyCalibrator()
|
||||
assert not cal.is_calibrated
|
||||
|
||||
def test_get_linear_calibration(self):
|
||||
cal = EnergyCalibrator()
|
||||
result = cal.get_linear_calibration(
|
||||
low_energy_keV=0.0,
|
||||
high_energy_keV=100.0,
|
||||
low_channel=0,
|
||||
high_channel=100,
|
||||
)
|
||||
assert cal.is_calibrated
|
||||
assert result.order == 1
|
||||
np.testing.assert_allclose(result(50), 50.0)
|
||||
|
||||
def test_apply_calibration_with_points(self):
|
||||
cal = EnergyCalibrator()
|
||||
cal.add_calibration_point(0, 0.0)
|
||||
cal.add_calibration_point(50, 50.0)
|
||||
cal.add_calibration_point(100, 100.0)
|
||||
|
||||
result = cal.fit_calibration(order=1)
|
||||
assert result.order == 1
|
||||
# Perfect fit on collinear points
|
||||
np.testing.assert_allclose(result(50), 50.0, rtol=1e-6)
|
||||
|
||||
def test_channel_to_energy_raises_without_calibration(self):
|
||||
cal = EnergyCalibrator()
|
||||
with pytest.raises(RuntimeError):
|
||||
cal.channel_to_energy(50)
|
||||
|
||||
def test_energy_to_channel_raises_without_calibration(self):
|
||||
cal = EnergyCalibrator()
|
||||
with pytest.raises(RuntimeError):
|
||||
cal.energy_to_channel(50.0)
|
||||
|
||||
def test_full_channel_mapping(self):
|
||||
"""Calibrate using full channel-to-energy mapping."""
|
||||
# Simulate: E = 0.5 + 0.1*ch + 0.001*ch²
|
||||
n_channels = 200
|
||||
channels = np.arange(n_channels, dtype=np.float64)
|
||||
energies = 0.5 + 0.1 * channels + 0.001 * channels**2
|
||||
|
||||
cal = EnergyCalibrator()
|
||||
result = cal.fit_calibration(order=2, channel_energies=energies)
|
||||
|
||||
assert result.order == 2
|
||||
np.testing.assert_allclose(result(channels), energies, rtol=1e-10)
|
||||
|
||||
def test_quadratic_calibration(self):
|
||||
cal = EnergyCalibrator()
|
||||
cal.add_calibration_point(0, 0.0)
|
||||
cal.add_calibration_point(50, 52.0)
|
||||
cal.add_calibration_point(100, 108.0) # slight quadratic deviation
|
||||
|
||||
result = cal.fit_calibration(order=2)
|
||||
assert result.order == 2
|
||||
@@ -0,0 +1,106 @@
|
||||
"""Tests for exp_data.corrections."""
|
||||
|
||||
import numpy as np
|
||||
|
||||
from exp_data.corrections import (
|
||||
DeadTimeCorrection,
|
||||
PileUpCorrection,
|
||||
BackgroundSubtractor,
|
||||
)
|
||||
|
||||
|
||||
class TestDeadTimeCorrection:
|
||||
"""Test dead-time correction."""
|
||||
|
||||
def test_non_paralyzable_live_time(self):
|
||||
dt = DeadTimeCorrection(dead_time=0.001, model="non_paralyzable")
|
||||
# At 100 cps: live = 1 - 100*0.001 = 0.9
|
||||
assert np.isclose(dt.live_time_fraction(100.0), 0.9)
|
||||
|
||||
def test_non_paralyzable_correct_rate(self):
|
||||
dt = DeadTimeCorrection(dead_time=0.01, model="non_paralyzable")
|
||||
# True rate = observed / (1 - obs*dt) = 90 / (1 - 90*0.01) = 90 / 0.1 = 900
|
||||
assert np.isclose(dt.correct_rate(90.0), 900.0, rtol=0.01)
|
||||
|
||||
def test_paralyzable_live_time(self):
|
||||
dt = DeadTimeCorrection(dead_time=0.001, model="paralyzable")
|
||||
# At low rate: exp(-0.1) ≈ 0.905
|
||||
frac = dt.live_time_fraction(100.0)
|
||||
assert np.isclose(frac, np.exp(-0.1), rtol=1e-5)
|
||||
|
||||
def test_correct_counts(self):
|
||||
dt = DeadTimeCorrection(dead_time=0.01, model="non_paralyzable")
|
||||
counts = np.array([100.0, 200.0, 300.0])
|
||||
live_time = 10.0
|
||||
corrected = dt.correct_counts(counts, live_time)
|
||||
# Total observed rate = 600/10 = 60 cps
|
||||
# True rate = 60 / (1 - 60*0.01) = 60 / 0.94 = 63.83 cps
|
||||
# Scale factor = 63.83/60 = 1.0638
|
||||
scale = 1.0 / (1 - 60.0 * 0.01)
|
||||
expected = counts * scale
|
||||
np.testing.assert_allclose(corrected, expected)
|
||||
|
||||
|
||||
class TestPileUpCorrection:
|
||||
"""Test pile-up correction."""
|
||||
|
||||
def test_pile_up_fraction(self):
|
||||
pu = PileUpCorrection(dead_time=0.001)
|
||||
# At 100 cps: (0.1)² = 0.01
|
||||
frac = pu.estimate_pile_up_fraction(100.0)
|
||||
assert np.isclose(frac, 0.01)
|
||||
|
||||
def test_correct_spectrum(self):
|
||||
pu = PileUpCorrection(dead_time=0.001)
|
||||
counts = np.array([100.0, 200.0, 300.0])
|
||||
live_time = 10.0
|
||||
corrected = pu.correct_spectrum(counts, live_time)
|
||||
# Pile-up correction adds a small fraction back
|
||||
assert np.all(corrected >= counts)
|
||||
|
||||
|
||||
class TestBackgroundSubtractor:
|
||||
"""Test background subtraction."""
|
||||
|
||||
def test_constant_background(self):
|
||||
energies = np.linspace(0, 100, 200)
|
||||
counts = 50 * np.exp(-0.5 * ((energies - 50) / 5) ** 2) + 10.0 # peak + bg
|
||||
|
||||
corrected = BackgroundSubtractor.constant_background(
|
||||
counts,
|
||||
energies,
|
||||
low_roi=(0.0, 10.0),
|
||||
high_roi=(90.0, 100.0),
|
||||
)
|
||||
|
||||
# Peak should still be visible, background removed
|
||||
assert np.max(corrected) > 0
|
||||
# Background regions should be near zero
|
||||
assert np.mean(corrected[energies <= 10]) < 5
|
||||
assert np.mean(corrected[energies >= 90]) < 5
|
||||
|
||||
def test_polynomial_background(self):
|
||||
energies = np.linspace(0, 100, 200)
|
||||
counts = 50 * np.exp(-0.5 * ((energies - 50) / 5) ** 2) + 10.0 + 0.05 * energies
|
||||
|
||||
corrected = BackgroundSubtractor.polynomial_background(
|
||||
counts,
|
||||
energies,
|
||||
order=1,
|
||||
low_roi=(0.0, 10.0),
|
||||
high_roi=(90.0, 100.0),
|
||||
)
|
||||
|
||||
assert np.max(corrected) > 0
|
||||
|
||||
def test_roi_average(self):
|
||||
counts = np.array([10.0, 12.0, 50.0, 10.0, 11.0])
|
||||
mask = np.array([True, True, False, True, True])
|
||||
avg = BackgroundSubtractor.roi_average(counts, mask)
|
||||
np.testing.assert_allclose(avg, 10.75)
|
||||
|
||||
def test_roi_average_empty(self):
|
||||
counts = np.array([10.0, 20.0])
|
||||
mask = np.array([False, False])
|
||||
avg = BackgroundSubtractor.roi_average(counts, mask)
|
||||
assert avg == 0.0
|
||||
@@ -0,0 +1,114 @@
|
||||
"""Tests for exp_data.fitters."""
|
||||
|
||||
import numpy as np
|
||||
|
||||
from exp_data.fitters import GaussianFitter, PeakFitter
|
||||
|
||||
|
||||
class TestGaussianFitter:
|
||||
"""Test Gaussian peak fitting."""
|
||||
|
||||
def _make_synthetic_peak(
|
||||
self,
|
||||
center=50.0,
|
||||
sigma=2.0,
|
||||
amplitude=100.0,
|
||||
background=10.0,
|
||||
n_points=500,
|
||||
energy_range=(0, 100),
|
||||
):
|
||||
"""Create synthetic Gaussian peak with noise."""
|
||||
energies = np.linspace(energy_range[0], energy_range[1], n_points)
|
||||
counts = (
|
||||
amplitude * np.exp(-0.5 * ((energies - center) / sigma) ** 2) + background
|
||||
)
|
||||
counts += np.random.RandomState(42).normal(
|
||||
0, np.sqrt(np.maximum(counts, 1)), n_points
|
||||
)
|
||||
return energies, counts
|
||||
|
||||
def test_fit_returns_valid_result(self):
|
||||
energies, counts = self._make_synthetic_peak()
|
||||
fitter = GaussianFitter(energies, counts)
|
||||
result = fitter.fit()
|
||||
|
||||
assert result.success
|
||||
assert result.chi2 > 0
|
||||
assert result.n_dof > 0
|
||||
assert result.covariance.shape == (4, 4)
|
||||
|
||||
def test_fit_recovers_parameters(self):
|
||||
center, sigma, amplitude, background = 50.0, 2.0, 100.0, 10.0
|
||||
energies, counts = self._make_synthetic_peak(
|
||||
center=center, sigma=sigma, amplitude=amplitude, background=background
|
||||
)
|
||||
fitter = GaussianFitter(energies, counts)
|
||||
result = fitter.fit()
|
||||
|
||||
np.testing.assert_allclose(result.center_keV, center, rtol=0.05)
|
||||
np.testing.assert_allclose(result.sigma_keV, sigma, rtol=0.1)
|
||||
np.testing.assert_allclose(result.amplitude, amplitude, rtol=0.15)
|
||||
|
||||
def test_resolution_property(self):
|
||||
energies, counts = self._make_synthetic_peak(sigma=2.0)
|
||||
fitter = GaussianFitter(energies, counts)
|
||||
result = fitter.fit()
|
||||
|
||||
expected_fwhm = 2.35482 * 2.0
|
||||
np.testing.assert_allclose(result.resolution(), expected_fwhm, rtol=0.05)
|
||||
|
||||
def test_resolution_fraction(self):
|
||||
energies, counts = self._make_synthetic_peak(center=50.0, sigma=2.0)
|
||||
fitter = GaussianFitter(energies, counts)
|
||||
result = fitter.fit()
|
||||
|
||||
expected_frac = (2.35482 * 2.0) / 50.0
|
||||
np.testing.assert_allclose(
|
||||
result.resolution_fraction(), expected_frac, rtol=0.05
|
||||
)
|
||||
|
||||
def test_auto_guess(self):
|
||||
energies, counts = self._make_synthetic_peak()
|
||||
fitter = GaussianFitter(energies, counts)
|
||||
result = fitter.fit()
|
||||
assert result.success
|
||||
|
||||
|
||||
class TestPeakFitter:
|
||||
"""Test multi-peak detection and fitting."""
|
||||
|
||||
def _make_multi_peak(self):
|
||||
"""Create spectrum with two peaks."""
|
||||
energies = np.linspace(0, 100, 1000)
|
||||
counts = np.zeros_like(energies)
|
||||
counts += 50 * np.exp(-0.5 * ((energies - 20) / 2) ** 2)
|
||||
counts += 80 * np.exp(-0.5 * ((energies - 70) / 3) ** 2)
|
||||
counts += 5 # constant background
|
||||
counts += np.random.RandomState(42).normal(
|
||||
0, np.sqrt(np.maximum(counts, 1)), len(energies)
|
||||
)
|
||||
return energies, counts
|
||||
|
||||
def test_find_peaks(self):
|
||||
energies, counts = self._make_multi_peak()
|
||||
fitter = PeakFitter(energies, counts)
|
||||
peaks = fitter.find_peaks(min_height=10.0, min_distance_keV=10.0)
|
||||
|
||||
assert len(peaks) >= 1 # At least one peak should be found
|
||||
|
||||
def test_fit_peaks_with_known_energies(self):
|
||||
energies, counts = self._make_multi_peak()
|
||||
fitter = PeakFitter(energies, counts)
|
||||
results = fitter.fit_peaks(known_energies=np.array([20.0, 70.0]))
|
||||
|
||||
assert len(results) >= 1
|
||||
for r in results:
|
||||
assert r.success
|
||||
assert r.n_dof > 0
|
||||
|
||||
def test_fit_peaks_auto_detect(self):
|
||||
energies, counts = self._make_multi_peak()
|
||||
fitter = PeakFitter(energies, counts)
|
||||
results = fitter.fit_peaks()
|
||||
|
||||
assert len(results) >= 1
|
||||
@@ -0,0 +1,141 @@
|
||||
"""Tests for exp_data.spectrum."""
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from exp_data.spectrum import ExpSpectrum
|
||||
|
||||
|
||||
class TestExpSpectrumBasic:
|
||||
"""Basic ExpSpectrum construction and properties."""
|
||||
|
||||
def test_construct_with_errors(self):
|
||||
spectrum = ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0, 2.0]),
|
||||
counts=np.array([10.0, 20.0, 30.0]),
|
||||
errors=np.array([3.16, 4.47, 5.48]),
|
||||
)
|
||||
assert spectrum.n_channels == 3
|
||||
np.testing.assert_array_equal(spectrum.energies, [0.0, 1.0, 2.0])
|
||||
|
||||
def test_automatic_poisson_errors(self):
|
||||
spectrum = ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0, 2.0]),
|
||||
counts=np.array([10.0, 20.0, 30.0]),
|
||||
)
|
||||
expected = np.sqrt([10.0, 20.0, 30.0])
|
||||
np.testing.assert_allclose(spectrum.errors, expected)
|
||||
|
||||
def test_zero_counts_get_error_one(self):
|
||||
spectrum = ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0]),
|
||||
counts=np.array([0.0, 100.0]),
|
||||
)
|
||||
assert spectrum.errors[0] == 1.0 # sqrt(1) for zero counts
|
||||
|
||||
def test_shape_mismatch_raises(self):
|
||||
with pytest.raises(AssertionError):
|
||||
ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0]),
|
||||
counts=np.array([10.0, 20.0, 30.0]),
|
||||
)
|
||||
|
||||
def test_negative_energies_raises(self):
|
||||
with pytest.raises(AssertionError):
|
||||
ExpSpectrum(
|
||||
energies=np.array([-1.0, 0.0]),
|
||||
counts=np.array([10.0, 20.0]),
|
||||
)
|
||||
|
||||
def test_negative_counts_raises(self):
|
||||
with pytest.raises(AssertionError):
|
||||
ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0]),
|
||||
counts=np.array([10.0, -5.0]),
|
||||
)
|
||||
|
||||
def test_total_counts(self):
|
||||
spectrum = ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0, 2.0]),
|
||||
counts=np.array([10.0, 20.0, 30.0]),
|
||||
)
|
||||
assert spectrum.total_counts == 60.0
|
||||
|
||||
def test_energy_range(self):
|
||||
spectrum = ExpSpectrum(
|
||||
energies=np.array([0.5, 1.5, 2.5]),
|
||||
counts=np.array([10.0, 20.0, 30.0]),
|
||||
)
|
||||
assert spectrum.energy_range == (0.5, 2.5)
|
||||
|
||||
def test_normalize(self):
|
||||
spectrum = ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0, 2.0]),
|
||||
counts=np.array([10.0, 20.0, 30.0]),
|
||||
)
|
||||
norm = spectrum.normalize()
|
||||
integral = np.trapezoid(norm, spectrum.energies)
|
||||
assert np.isclose(integral, 1.0, rtol=1e-10)
|
||||
|
||||
|
||||
class TestExpSpectrumDeadTime:
|
||||
"""Dead-time correction on ExpSpectrum."""
|
||||
|
||||
def test_no_correction_when_dead_time_zero(self):
|
||||
spectrum = ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0]),
|
||||
counts=np.array([10.0, 20.0]),
|
||||
)
|
||||
corrected = spectrum.apply_dead_time_correction()
|
||||
np.testing.assert_array_equal(corrected.counts, spectrum.counts)
|
||||
|
||||
def test_correction_applied(self):
|
||||
spectrum = ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0]),
|
||||
counts=np.array([10.0, 20.0]),
|
||||
dead_time=0.1,
|
||||
)
|
||||
corrected = spectrum.apply_dead_time_correction()
|
||||
expected = spectrum.counts / 0.9
|
||||
np.testing.assert_allclose(corrected.counts, expected, rtol=1e-10)
|
||||
assert corrected.dead_time == 0.0
|
||||
|
||||
def test_live_time_corrected_rate(self):
|
||||
spectrum = ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0]),
|
||||
counts=np.array([10.0, 20.0]),
|
||||
dead_time=0.1,
|
||||
live_time=10.0,
|
||||
)
|
||||
rate = spectrum.live_time_corrected_rate
|
||||
assert rate is not None
|
||||
# Total counts / (live_time * (1 - dead_time))
|
||||
expected = 30.0 / (10.0 * 0.9)
|
||||
assert np.isclose(rate, expected)
|
||||
|
||||
|
||||
class TestExpSpectrumSerialization:
|
||||
"""Dict round-trip for ExpSpectrum."""
|
||||
|
||||
def test_to_from_dict(self):
|
||||
spectrum = ExpSpectrum(
|
||||
energies=np.array([0.0, 1.0, 2.0]),
|
||||
counts=np.array([10.0, 20.0, 30.0]),
|
||||
dead_time=0.05,
|
||||
live_time=100.0,
|
||||
source="Tc99",
|
||||
run_id="run_042",
|
||||
date="2025-01-15",
|
||||
metadata={"detector": "MMC", "temperature": "25 mK"},
|
||||
)
|
||||
|
||||
recovered = ExpSpectrum.from_dict(spectrum.to_dict())
|
||||
|
||||
np.testing.assert_allclose(recovered.energies, spectrum.energies)
|
||||
np.testing.assert_allclose(recovered.counts, spectrum.counts)
|
||||
assert recovered.dead_time == spectrum.dead_time
|
||||
assert recovered.live_time == spectrum.live_time
|
||||
assert recovered.source == spectrum.source
|
||||
assert recovered.run_id == spectrum.run_id
|
||||
assert recovered.date == spectrum.date
|
||||
assert recovered.metadata == spectrum.metadata
|
||||
@@ -0,0 +1,161 @@
|
||||
"""Tests for fitter.extractor."""
|
||||
|
||||
import numpy as np
|
||||
|
||||
from fitter.extractor import (
|
||||
CWExtractor,
|
||||
GVAExtractor,
|
||||
CWExtractionResult,
|
||||
GVAExtractionResult,
|
||||
)
|
||||
|
||||
|
||||
class TestCWExtractor:
|
||||
"""Test CWExtractor."""
|
||||
|
||||
def test_extract(self):
|
||||
from beta_spectrum import SpectrumConfig
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
)
|
||||
extractor = CWExtractor(config)
|
||||
|
||||
energies_keV = np.linspace(0.01, 0.29, 50)
|
||||
measured = np.linspace(10.0, 0.1, 50)
|
||||
errors = np.sqrt(measured)
|
||||
model = np.linspace(9.5, 0.09, 50)
|
||||
|
||||
result = extractor.extract(
|
||||
measured_counts=measured,
|
||||
measured_errors=errors,
|
||||
energies_keV=energies_keV,
|
||||
model_counts=model,
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
assert isinstance(result, CWExtractionResult)
|
||||
assert len(result.cw_values) == 50
|
||||
assert len(result.energies_keV) == 50
|
||||
assert np.isclose(result.endpoint_keV, 294.0)
|
||||
assert np.all(np.isfinite(result.cw_values))
|
||||
|
||||
def test_extract_with_zero_model(self):
|
||||
from beta_spectrum import SpectrumConfig
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
)
|
||||
extractor = CWExtractor(config)
|
||||
|
||||
energies_keV = np.linspace(0.01, 0.29, 50)
|
||||
measured = np.ones(50) * 5.0
|
||||
errors = np.ones(50) * 2.0
|
||||
model = np.zeros(50) # All zeros
|
||||
|
||||
result = extractor.extract(
|
||||
measured_counts=measured,
|
||||
measured_errors=errors,
|
||||
energies_keV=energies_keV,
|
||||
model_counts=model,
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
# C(W) should be 1.0 where model is zero (default)
|
||||
assert np.all(result.cw_values == 1.0)
|
||||
|
||||
|
||||
class TestGVAExtractor:
|
||||
"""Test GVAExtractor."""
|
||||
|
||||
def test_extract_from_curve_basic(self):
|
||||
energies_W = np.linspace(1.1, 1.57, 100)
|
||||
endpoint_W = 1.575 # ~294 keV
|
||||
cw_values = np.ones_like(energies_W) # C(W) = 1 for allowed
|
||||
|
||||
result = GVAExtractor.extract_from_curve(
|
||||
cw_values=cw_values,
|
||||
energies_W=energies_W,
|
||||
endpoint_W=endpoint_W,
|
||||
)
|
||||
|
||||
assert isinstance(result, GVAExtractionResult)
|
||||
assert np.isclose(result.gv_eff, 1.0)
|
||||
assert result.chi2_per_dof >= 0
|
||||
|
||||
def test_extract_from_curve_with_slope(self):
|
||||
energies_W = np.linspace(1.1, 1.57, 100)
|
||||
endpoint_W = 1.575
|
||||
# C(W) with slight slope
|
||||
cw_values = 1.0 + 0.1 * (energies_W - endpoint_W)
|
||||
|
||||
result = GVAExtractor.extract_from_curve(
|
||||
cw_values=cw_values,
|
||||
energies_W=energies_W,
|
||||
endpoint_W=endpoint_W,
|
||||
)
|
||||
|
||||
assert isinstance(result, GVAExtractionResult)
|
||||
# g_A should be non-zero due to slope
|
||||
assert abs(result.ga_eff) > 0
|
||||
|
||||
def test_extract_from_curve_insufficient_data(self):
|
||||
energies_W = np.array([1.0, 1.5])
|
||||
endpoint_W = 1.575
|
||||
cw_values = np.array([1.0, 1.0])
|
||||
|
||||
result = GVAExtractor.extract_from_curve(
|
||||
cw_values=cw_values,
|
||||
energies_W=energies_W,
|
||||
endpoint_W=endpoint_W,
|
||||
)
|
||||
|
||||
assert isinstance(result, GVAExtractionResult)
|
||||
assert result.chi2_per_dof == 0.0
|
||||
|
||||
def test_gva_summary(self):
|
||||
result = GVAExtractionResult(
|
||||
gv_eff=1.0,
|
||||
gv_error=0.01,
|
||||
ga_eff=0.5,
|
||||
ga_error=0.02,
|
||||
correlation=-0.3,
|
||||
chi2_per_dof=1.2,
|
||||
)
|
||||
summary = result.summary()
|
||||
assert "g_V" in summary
|
||||
assert "g_A" in summary
|
||||
assert "1.0" in summary
|
||||
assert "0.5" in summary
|
||||
assert "1.2" in summary
|
||||
|
||||
|
||||
class TestCWExtractionResult:
|
||||
"""Test CWExtractionResult methods."""
|
||||
|
||||
def test_kurie_plot(self):
|
||||
energies_W = np.linspace(1.1, 1.57, 100)
|
||||
endpoint_W = 1.575
|
||||
|
||||
result = CWExtractionResult(
|
||||
energies_W=energies_W,
|
||||
energies_keV=energies_W * 510.998950 - 510.998950,
|
||||
cw_values=np.ones_like(energies_W),
|
||||
cw_errors=np.zeros_like(energies_W),
|
||||
flux=np.ones_like(energies_W) * 10.0,
|
||||
flux_errors=np.ones_like(energies_W) * 3.0,
|
||||
endpoint_W=endpoint_W,
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
E_keV, kurie = result.kurie_plot()
|
||||
assert len(E_keV) > 0
|
||||
assert len(kurie) > 0
|
||||
assert np.all(np.isfinite(kurie))
|
||||
assert np.all(kurie >= 0)
|
||||
@@ -0,0 +1,176 @@
|
||||
"""Tests for fitter.fit_engine."""
|
||||
|
||||
import numpy as np
|
||||
|
||||
from fitter.fit_engine import SpectrumFitter, FitConfig
|
||||
from fitter.model import SpectrumModel
|
||||
|
||||
|
||||
class TestFitConfig:
|
||||
"""Test FitConfig defaults."""
|
||||
|
||||
def test_defaults(self):
|
||||
config = FitConfig()
|
||||
assert config.method == "trf"
|
||||
assert config.max_nfev == 1000
|
||||
assert config.loss == "linear"
|
||||
|
||||
|
||||
class TestSpectrumFitter:
|
||||
"""Test SpectrumFitter."""
|
||||
|
||||
def _make_fitter(self):
|
||||
"""Create a simple fitter with synthetic data."""
|
||||
|
||||
# Create a simple linear model
|
||||
def linear_model(energies_keV, norm, bg):
|
||||
return norm * (1.0 - energies_keV / 100.0) + bg
|
||||
|
||||
# Synthetic data
|
||||
energies = np.linspace(0.1, 99.9, 200)
|
||||
counts = 10.0 * (1.0 - energies / 100.0) + 2.0
|
||||
counts += np.random.RandomState(42).normal(0, 0.5, len(energies))
|
||||
|
||||
# Mock model
|
||||
class MockModel:
|
||||
def evaluate(self, energies_keV, norm, bg):
|
||||
return linear_model(energies_keV, norm, bg)
|
||||
|
||||
model = MockModel()
|
||||
fitter = SpectrumFitter(
|
||||
model, energies, counts, exp_errors=np.ones_like(counts) * 0.5
|
||||
)
|
||||
return fitter
|
||||
|
||||
def test_fit_basic(self):
|
||||
fitter = self._make_fitter()
|
||||
result = fitter.fit(x0=[10.0, 2.0], param_names=["norm", "bg"])
|
||||
|
||||
assert result.success
|
||||
assert result.n_points == 200
|
||||
assert result.n_free == 2
|
||||
assert result.chi2 > 0
|
||||
assert result.chi2_per_dof > 0
|
||||
assert len(result.parameters) == 2
|
||||
|
||||
def test_fit_with_bounds(self):
|
||||
fitter = self._make_fitter()
|
||||
result = fitter.fit(
|
||||
x0=[10.0, 2.0],
|
||||
bounds=([1.0, 0.0], [20.0, 10.0]),
|
||||
)
|
||||
assert result.success
|
||||
|
||||
def test_fit_returns_residuals(self):
|
||||
fitter = self._make_fitter()
|
||||
result = fitter.fit(x0=[10.0, 2.0])
|
||||
|
||||
assert len(result.residuals) == 200
|
||||
assert len(result.model_values) == 200
|
||||
|
||||
def test_parameters_property(self):
|
||||
fitter = self._make_fitter()
|
||||
fitter.fit(x0=[10.0, 2.0])
|
||||
params = fitter.parameters
|
||||
assert len(params) == 2
|
||||
|
||||
def test_covariance_property(self):
|
||||
fitter = self._make_fitter()
|
||||
fitter.fit(x0=[10.0, 2.0])
|
||||
cov = fitter.covariance
|
||||
assert cov.shape == (2, 2)
|
||||
|
||||
def test_chi2_per_dof(self):
|
||||
fitter = self._make_fitter()
|
||||
result = fitter.fit(x0=[10.0, 2.0])
|
||||
assert result.chi2_per_dof == result.chi2 / result.n_free
|
||||
|
||||
def test_p_value(self):
|
||||
fitter = self._make_fitter()
|
||||
result = fitter.fit(x0=[10.0, 2.0])
|
||||
assert 0.0 <= result.p_value <= 1.0
|
||||
|
||||
def test_parameters_with_errors(self):
|
||||
fitter = self._make_fitter()
|
||||
result = fitter.fit(x0=[10.0, 2.0])
|
||||
params = result.parameters_with_errors
|
||||
assert len(params) == 2
|
||||
for name, (val, err) in params.items():
|
||||
assert isinstance(val, float)
|
||||
assert isinstance(err, float)
|
||||
assert err >= 0
|
||||
|
||||
def test_correlation_matrix(self):
|
||||
fitter = self._make_fitter()
|
||||
result = fitter.fit(x0=[10.0, 2.0])
|
||||
corr = result.correlation_matrix
|
||||
assert corr.shape == (2, 2)
|
||||
# Diagonal should be 1
|
||||
np.testing.assert_allclose(np.diag(corr), 1.0)
|
||||
|
||||
def test_summary(self):
|
||||
fitter = self._make_fitter()
|
||||
result = fitter.fit(x0=[10.0, 2.0], param_names=["norm", "bg"])
|
||||
summary = result.summary(param_names=["norm", "bg"])
|
||||
assert "norm" in summary
|
||||
assert "bg" in summary
|
||||
assert "Fit Results" in summary
|
||||
|
||||
def test_profile_likelihood(self):
|
||||
fitter = self._make_fitter()
|
||||
result = fitter.fit(x0=[10.0, 2.0])
|
||||
param_vals, chi2_vals, norm_vals = result.profile_likelihood(0, n_points=10)
|
||||
assert len(param_vals) == 10
|
||||
assert len(chi2_vals) == 10
|
||||
assert len(norm_vals) == 10
|
||||
assert np.min(norm_vals) >= 0 # Normalized χ² should be >= 0
|
||||
|
||||
def test_confidence_interval(self):
|
||||
fitter = self._make_fitter()
|
||||
result = fitter.fit(x0=[10.0, 2.0])
|
||||
ci = result.confidence_interval(0, confidence=0.95)
|
||||
assert len(ci) == 2
|
||||
assert ci[0] <= ci[1]
|
||||
|
||||
|
||||
class TestSpectrumFitterWithRealModel:
|
||||
"""Test SpectrumFitter with a real SpectrumModel."""
|
||||
|
||||
def test_fit_with_real_model(self):
|
||||
from beta_spectrum import BetaSpectrum, SpectrumConfig
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
detector_model="gaussian",
|
||||
detector_sigma_a_keV=0.5,
|
||||
detector_n_channels=200,
|
||||
)
|
||||
spectrum = BetaSpectrum.from_config(config)
|
||||
detector = BetaSpectrum.create_detector_from_config(config)
|
||||
model = SpectrumModel(spectrum, detector)
|
||||
|
||||
# Generate synthetic data from model
|
||||
energies_keV = detector.channel_energies
|
||||
true_norm = 5.0
|
||||
true_bg = 0.1
|
||||
counts = model.evaluate(energies_keV, true_norm, true_bg)
|
||||
counts = np.maximum(counts, 0) # Clip negatives
|
||||
counts += np.random.RandomState(42).normal(
|
||||
0, np.sqrt(np.maximum(counts, 1)), len(energies_keV)
|
||||
)
|
||||
counts = np.maximum(counts, 0)
|
||||
|
||||
fitter = SpectrumFitter(
|
||||
model,
|
||||
energies_keV,
|
||||
counts,
|
||||
exp_errors=np.sqrt(np.maximum(counts, 1)),
|
||||
)
|
||||
result = fitter.fit(x0=[1.0, 0.0], param_names=["norm", "bg"])
|
||||
|
||||
assert result.success
|
||||
# Fit should complete without error
|
||||
assert len(result.parameters) == 2
|
||||
@@ -0,0 +1,119 @@
|
||||
"""Tests for fitter.model."""
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from beta_spectrum import BetaSpectrum, SpectrumConfig
|
||||
from fitter.model import SpectrumModel
|
||||
|
||||
|
||||
class TestSpectrumModel:
|
||||
"""Test SpectrumModel."""
|
||||
|
||||
def _make_model(self, with_detector: bool = False):
|
||||
"""Create a test SpectrumModel."""
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
)
|
||||
spectrum = BetaSpectrum.from_config(config)
|
||||
|
||||
detector = None
|
||||
if with_detector:
|
||||
detector = BetaSpectrum.create_detector_from_config(
|
||||
SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
detector_model="gaussian",
|
||||
detector_sigma_a_keV=1.0,
|
||||
detector_n_channels=100,
|
||||
)
|
||||
)
|
||||
|
||||
model = SpectrumModel(spectrum, detector)
|
||||
return model, config
|
||||
|
||||
def test_evaluate_unconvolved(self):
|
||||
model, config = self._make_model()
|
||||
energies_keV = np.linspace(0.01, 0.29, 50)
|
||||
values = model.evaluate_unconvolved(energies_keV)
|
||||
assert len(values) == 50
|
||||
assert np.all(np.isfinite(values))
|
||||
assert np.all(values >= 0)
|
||||
|
||||
def test_evaluate_no_detector(self):
|
||||
model, config = self._make_model()
|
||||
energies_keV = np.linspace(0.01, 0.29, 50)
|
||||
values = model.evaluate(energies_keV, normalization=1.0, background=0.0)
|
||||
assert len(values) == 50
|
||||
assert np.all(np.isfinite(values))
|
||||
|
||||
def test_evaluate_with_normalization(self):
|
||||
model, config = self._make_model()
|
||||
energies_keV = np.linspace(0.01, 0.29, 50)
|
||||
values1 = model.evaluate(energies_keV, normalization=1.0, background=0.0)
|
||||
values2 = model.evaluate(energies_keV, normalization=2.0, background=0.0)
|
||||
np.testing.assert_allclose(values2, values1 * 2.0)
|
||||
|
||||
def test_evaluate_with_background(self):
|
||||
model, config = self._make_model()
|
||||
energies_keV = np.linspace(0.01, 0.29, 50)
|
||||
values1 = model.evaluate(energies_keV, normalization=1.0, background=0.0)
|
||||
values2 = model.evaluate(energies_keV, normalization=1.0, background=5.0)
|
||||
np.testing.assert_allclose(values2, values1 + 5.0)
|
||||
|
||||
def test_convolve_requires_detector(self):
|
||||
model, _ = self._make_model()
|
||||
energies_keV = np.linspace(0.01, 0.29, 50)
|
||||
with pytest.raises(RuntimeError):
|
||||
model.convolve(energies_keV)
|
||||
|
||||
def test_convolve_with_detector(self):
|
||||
model, _ = self._make_model(with_detector=True)
|
||||
energies_keV = model.detector_response.channel_energies
|
||||
values = model.convolve(energies_keV)
|
||||
assert len(values) == model.detector_response.n_channels
|
||||
assert np.all(np.isfinite(values))
|
||||
# Convolution can produce small negative values due to numerical effects
|
||||
# — check that most values are non-negative
|
||||
assert np.sum(values >= 0) > len(values) * 0.9
|
||||
|
||||
def test_default_params(self):
|
||||
from exp_data.spectrum import ExpSpectrum
|
||||
|
||||
energies = np.linspace(0, 100, 200)
|
||||
counts = 50 * np.exp(-0.5 * ((energies - 50) / 5) ** 2) + 10
|
||||
exp_spec = ExpSpectrum(energies=energies, counts=counts)
|
||||
|
||||
norm, bg = SpectrumModel.default_params(exp_spec)
|
||||
assert norm > 0
|
||||
assert bg >= 0
|
||||
|
||||
def test_default_bounds(self):
|
||||
lower, upper = SpectrumModel.default_bounds(100)
|
||||
assert len(lower) == 2
|
||||
assert len(upper) == 2
|
||||
assert lower[0] < upper[0]
|
||||
assert lower[1] < upper[1]
|
||||
|
||||
|
||||
class TestSpectrumModelKeVToW:
|
||||
"""Test energy unit conversion in SpectrumModel."""
|
||||
|
||||
def test_keV_to_W(self):
|
||||
model = SpectrumModel(None) # type: ignore
|
||||
# 0 keV → W = 1.0
|
||||
W = model._keV_to_W(np.array([0.0]))
|
||||
assert np.isclose(W[0], 1.0)
|
||||
|
||||
# 511 keV → W = 2.0
|
||||
W = model._keV_to_W(np.array([510.998950]))
|
||||
assert np.isclose(W[0], 2.0)
|
||||
|
||||
# 294 keV → W ≈ 1.575
|
||||
W = model._keV_to_W(np.array([294.0]))
|
||||
assert np.isclose(W[0], 294.0 / 510.998950 + 1.0)
|
||||
@@ -0,0 +1,100 @@
|
||||
"""Tests for fitter.result."""
|
||||
|
||||
import numpy as np
|
||||
|
||||
from fitter.result import AnalysisResult, FitSummary
|
||||
|
||||
|
||||
class MockFitResult:
|
||||
"""Mock FitResult for testing AnalysisResult."""
|
||||
|
||||
def __init__(self):
|
||||
self.success = True
|
||||
self.chi2 = 12.5
|
||||
self.n_points = 100
|
||||
self.n_free = 3
|
||||
self.parameters = np.array([5.0, 0.1, 2.0])
|
||||
self.covariance = np.array(
|
||||
[
|
||||
[0.01, 0.0, 0.0],
|
||||
[0.0, 0.001, 0.0],
|
||||
[0.0, 0.0, 0.04],
|
||||
]
|
||||
)
|
||||
self.model_values = np.ones(100)
|
||||
self.residuals = np.zeros(100)
|
||||
self.message = "Optimization terminated successfully."
|
||||
|
||||
@property
|
||||
def chi2_per_dof(self):
|
||||
return self.chi2 / (self.n_points - self.n_free)
|
||||
|
||||
@property
|
||||
def p_value(self):
|
||||
from scipy.stats import chi2
|
||||
|
||||
return 1.0 - chi2.cdf(self.chi2, self.n_points - self.n_free)
|
||||
|
||||
@property
|
||||
def parameters_with_errors(self):
|
||||
errors = np.sqrt(np.diag(self.covariance))
|
||||
return {
|
||||
f"p{i}": (float(self.parameters[i]), float(errors[i]))
|
||||
for i in range(len(self.parameters))
|
||||
}
|
||||
|
||||
|
||||
class TestFitSummary:
|
||||
"""Test FitSummary."""
|
||||
|
||||
def test_str(self):
|
||||
summary = FitSummary(
|
||||
nuclide="Tc99",
|
||||
converged=True,
|
||||
chi2_per_dof=1.234,
|
||||
p_value=0.0567,
|
||||
parameters={"norm": (5.0, 0.1), "bg": (0.1, 0.01)},
|
||||
)
|
||||
text = str(summary)
|
||||
assert "Tc99" in text
|
||||
assert "1.234" in text
|
||||
assert "5.0" in text
|
||||
assert "Converged: True" in text
|
||||
|
||||
|
||||
class TestAnalysisResult:
|
||||
"""Test AnalysisResult."""
|
||||
|
||||
def test_summary(self):
|
||||
mock = MockFitResult()
|
||||
result = AnalysisResult(
|
||||
fit_result=mock,
|
||||
nuclide="Tc99",
|
||||
metadata={"param_names": ["norm", "bg", "endpoint"]},
|
||||
)
|
||||
summary = result.summary
|
||||
assert summary.nuclide == "Tc99"
|
||||
assert summary.converged is True
|
||||
assert summary.chi2_per_dof > 0
|
||||
|
||||
def test_to_dict(self):
|
||||
mock = MockFitResult()
|
||||
result = AnalysisResult(
|
||||
fit_result=mock,
|
||||
nuclide="Tc99",
|
||||
gv_eff=0.97,
|
||||
ga_eff=0.57,
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
d = result.to_dict()
|
||||
assert d["nuclide"] == "Tc99"
|
||||
assert d["gv_eff"] == 0.97
|
||||
assert d["ga_eff"] == 0.57
|
||||
assert d["endpoint_keV"] == 294.0
|
||||
assert d["converged"] is True
|
||||
|
||||
def test_empty_result(self):
|
||||
mock = MockFitResult()
|
||||
result = AnalysisResult(fit_result=mock)
|
||||
summary = result.summary
|
||||
assert summary.nuclide == ""
|
||||
@@ -1,675 +0,0 @@
|
||||
# tests/test_cw_extractor.py — Tests for C(W) shape factor extraction
|
||||
|
||||
"""
|
||||
Why test this?
|
||||
--------------
|
||||
The C(W) extractor is the core scientific tool for extracting the shape factor
|
||||
from measured beta spectra. It must correctly handle:
|
||||
1. Energy unit conversions (keV ↔ m_e units)
|
||||
2. Theoretical factor computation (Fermi, screening, radiative, etc.)
|
||||
3. Normalization of extracted C(W)
|
||||
4. Kurie plot endpoint determination
|
||||
5. g_V, g_A extraction from C(W) data
|
||||
|
||||
We test with synthetic data where the true answer is known, so we can verify
|
||||
the extraction recovers the correct values.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from beta_spectrum import (
|
||||
BetaSpectrum,
|
||||
CWExtractor,
|
||||
CWExtractionResult,
|
||||
GVAExtractionResult,
|
||||
SpectrumConfig,
|
||||
T_to_W,
|
||||
W_to_T,
|
||||
)
|
||||
|
||||
|
||||
class TestCWExtractionBasic:
|
||||
"""Basic C(W) extraction tests."""
|
||||
|
||||
def _make_synthetic_data(
|
||||
self,
|
||||
endpoint_keV: float = 294.0,
|
||||
n_points: int = 100,
|
||||
seed: int = 42,
|
||||
):
|
||||
"""Create synthetic beta spectrum with C(W) = 1 (allowed decay)."""
|
||||
np.random.seed(seed)
|
||||
|
||||
# Energy grid
|
||||
energies_keV = np.linspace(5, endpoint_keV * 0.98, n_points)
|
||||
|
||||
# Build theoretical spectrum with C(W) = 1
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=endpoint_keV / 1000.0,
|
||||
transition_type="F1",
|
||||
use_screening=True,
|
||||
use_exchange=True,
|
||||
use_finite_size=True,
|
||||
use_radiative=True,
|
||||
)
|
||||
spectrum = BetaSpectrum.from_config(config)
|
||||
W, E_MeV = spectrum.get_energy_grid(config)
|
||||
|
||||
# Interpolate to our energy grid
|
||||
from scipy.interpolate import interp1d
|
||||
|
||||
flux_interp = interp1d(
|
||||
E_MeV * 1000, spectrum(W), kind="linear", fill_value=0.0, bounds_error=False
|
||||
)
|
||||
flux_theory = flux_interp(energies_keV)
|
||||
|
||||
# Add noise (Poisson-like)
|
||||
noise_std = np.sqrt(np.maximum(flux_theory, 1e-10))
|
||||
flux = flux_theory + np.random.normal(0, noise_std, len(flux_theory))
|
||||
flux = np.maximum(flux, 0.0) # Flux cannot be negative
|
||||
flux_errors = np.sqrt(np.maximum(flux, 1e-10))
|
||||
|
||||
return energies_keV, flux, flux_errors, config
|
||||
|
||||
def test_CW_extraction_allowed_decay(self):
|
||||
"""For allowed decay (C(W)=1), extracted C(W) should be ~1."""
|
||||
energies_keV, flux, flux_errors, config = self._make_synthetic_data()
|
||||
|
||||
extractor = CWExtractor(config, flux, energies_keV, flux_errors)
|
||||
result = extractor.extract_CW(endpoint_keV=294.0)
|
||||
|
||||
# C(W) should be close to 1 in mid-energy region
|
||||
mid_mask = (energies_keV > 50) & (energies_keV < 200)
|
||||
if np.any(mid_mask):
|
||||
avg_cw = np.mean(result.cw_values[mid_mask])
|
||||
assert (
|
||||
abs(avg_cw - 1.0) < 0.3
|
||||
), f"Allowed decay C(W) should be ~1, got {avg_cw}"
|
||||
|
||||
def test_CW_extraction_preserves_endpoint(self):
|
||||
"""Extracted C(W) should report correct endpoint."""
|
||||
energies_keV, flux, flux_errors, config = self._make_synthetic_data()
|
||||
|
||||
extractor = CWExtractor(config, flux, energies_keV, flux_errors)
|
||||
result = extractor.extract_CW(endpoint_keV=294.0)
|
||||
|
||||
assert (
|
||||
abs(result.endpoint_keV - 294.0) < 1.0
|
||||
), f"Endpoint should be ~294 keV, got {result.endpoint_keV}"
|
||||
|
||||
def test_CW_extraction_errors_positive(self):
|
||||
"""C(W) uncertainties should be positive."""
|
||||
energies_keV, flux, flux_errors, config = self._make_synthetic_data()
|
||||
|
||||
extractor = CWExtractor(config, flux, energies_keV, flux_errors)
|
||||
result = extractor.extract_CW(endpoint_keV=294.0)
|
||||
|
||||
assert np.all(result.cw_errors > 0), "All errors should be positive"
|
||||
|
||||
def test_CW_extraction_clipped(self):
|
||||
"""C(W) values should be clipped to [0.01, 10.0]."""
|
||||
energies_keV, flux, flux_errors, config = self._make_synthetic_data()
|
||||
|
||||
extractor = CWExtractor(config, flux, energies_keV, flux_errors)
|
||||
result = extractor.extract_CW(endpoint_keV=294.0)
|
||||
|
||||
assert np.all(result.cw_values >= 0.01), "C(W) should be >= 0.01"
|
||||
assert np.all(result.cw_values <= 10.0), "C(W) should be <= 10.0"
|
||||
|
||||
|
||||
class TestCWParametrization:
|
||||
"""Test C(W) parametrization fitting."""
|
||||
|
||||
def _make_CW_data(
|
||||
self,
|
||||
endpoint_keV: float = 294.0,
|
||||
n_points: int = 60,
|
||||
seed: int = 42,
|
||||
):
|
||||
"""Create synthetic C(W) data with known parametrization."""
|
||||
np.random.seed(seed)
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=endpoint_keV / 1000.0,
|
||||
transition_type="F1",
|
||||
use_screening=True,
|
||||
use_exchange=True,
|
||||
use_finite_size=True,
|
||||
use_radiative=True,
|
||||
)
|
||||
|
||||
# Generate C(W) values with known form
|
||||
W = np.linspace(T_to_W(0.05), T_to_W(endpoint_keV * 0.95), n_points)
|
||||
E_keV = np.array([W_to_T(w) for w in W])
|
||||
|
||||
# Linear C(W): C(W) = 1.0 + 0.5*(W-1)
|
||||
true_a0, true_a1 = 1.0, 0.5
|
||||
cw_true = true_a0 + true_a1 * (W - 1.0)
|
||||
cw_noise = cw_true + np.random.normal(0, 0.05, len(W))
|
||||
|
||||
# C(W) = flux / (C0 * phase_space * Fermi * corrections * (W0-W)^2)
|
||||
# So: flux = C(W) * C0 * phase_space * Fermi * corrections * (W0-W)^2
|
||||
|
||||
p = np.sqrt(np.maximum(W**2 - 1.0, 0.0))
|
||||
from beta_spectrum.components.fermi import FermiFunction
|
||||
|
||||
F = FermiFunction(Z=config.Z_parent, A=config.A_number)(W)
|
||||
|
||||
C0 = 1.0 # Normalization
|
||||
endpoint_factor = (T_to_W(endpoint_keV / 1000.0) - W) ** 2
|
||||
|
||||
flux_synthetic = cw_noise * C0 * p * W * F * endpoint_factor
|
||||
flux_synthetic = np.maximum(flux_synthetic, 0)
|
||||
flux_errors = np.sqrt(np.maximum(flux_synthetic, 1.0))
|
||||
|
||||
return E_keV, flux_synthetic, flux_errors, config, true_a0, true_a1
|
||||
|
||||
def test_constant_parametrization(self):
|
||||
"""Fit C(W) = a₀ with known constant value."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(10, 250, 50)
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
cw_true = np.full_like(W, 1.0)
|
||||
cw = cw_true + np.random.normal(0, 0.05, len(W))
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
transition_type="F1",
|
||||
use_screening=True,
|
||||
use_exchange=True,
|
||||
use_finite_size=True,
|
||||
use_radiative=True,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(
|
||||
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.05
|
||||
)
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.05,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.05,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
result = extractor.fit_CW_parametrization(cw_result, parametrization="constant")
|
||||
|
||||
assert (
|
||||
abs(result.parameters[0] - 1.0) < 0.1
|
||||
), f"Constant C(W) should be ~1, got {result.parameters[0]}"
|
||||
|
||||
def test_linear_parametrization(self):
|
||||
"""Fit C(W) = a₀ + a₁·(W-1) with known coefficients."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(10, 250, 50)
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
true_a0, true_a1 = 1.0, 0.5
|
||||
cw = true_a0 + true_a1 * (W - 1.0) + np.random.normal(0, 0.05, len(W))
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
transition_type="F1",
|
||||
use_screening=True,
|
||||
use_exchange=True,
|
||||
use_finite_size=True,
|
||||
use_radiative=True,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(
|
||||
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.05
|
||||
)
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.05,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.05,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
result = extractor.fit_CW_parametrization(cw_result, parametrization="linear")
|
||||
|
||||
assert (
|
||||
abs(result.parameters[0] - true_a0) < 0.15
|
||||
), f"a₀: expected {true_a0}, got {result.parameters[0]}"
|
||||
assert (
|
||||
abs(result.parameters[1] - true_a1) < 0.5
|
||||
), f"a₁: expected {true_a1}, got {result.parameters[1]}"
|
||||
|
||||
def test_quadratic_parametrization(self):
|
||||
"""Fit C(W) = a₀ + a₁·(W-1) + a₂·(W-1)² with known coefficients."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(10, 250, 60)
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
true_a0, true_a1, true_a2 = 1.0, 0.3, -0.1
|
||||
cw = true_a0 + true_a1 * (W - 1.0) + true_a2 * (W - 1.0) ** 2
|
||||
cw += np.random.normal(0, 0.03, len(W))
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
transition_type="F1",
|
||||
use_screening=True,
|
||||
use_exchange=True,
|
||||
use_finite_size=True,
|
||||
use_radiative=True,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(
|
||||
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.03
|
||||
)
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.03,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.03,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
result = extractor.fit_CW_parametrization(
|
||||
cw_result, parametrization="quadratic"
|
||||
)
|
||||
|
||||
assert abs(result.parameters[0] - true_a0) < 0.15
|
||||
assert abs(result.parameters[1] - true_a1) < 0.5
|
||||
assert abs(result.parameters[2] - true_a2) < 0.5
|
||||
|
||||
def test_invalid_parametrization(self):
|
||||
"""Invalid parametrization should raise ValueError."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(10, 250, 30)
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
cw = np.ones_like(W)
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(config, cw, energies_keV)
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.05,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.05,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
with pytest.raises(ValueError, match="Unknown parametrization"):
|
||||
extractor.fit_CW_parametrization(cw_result, parametrization="invalid")
|
||||
|
||||
|
||||
class TestGVAExtraction:
|
||||
"""Test g_V, g_A extraction from C(W) data."""
|
||||
|
||||
def test_gV_gA_pure_fermi(self):
|
||||
"""For pure Fermi transition, g_A ≈ 0."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(10, 250, 50)
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
|
||||
# Pure Fermi: C(W) = g_V² (constant)
|
||||
true_gV = 0.5
|
||||
cw = true_gV**2 + np.random.normal(0, 0.02, len(W))
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
transition_type="F1",
|
||||
use_screening=True,
|
||||
use_exchange=True,
|
||||
use_finite_size=True,
|
||||
use_radiative=True,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(
|
||||
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.02
|
||||
)
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.02,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.02,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
result = extractor.fit_gV_gA(cw_result, M_F=1.0, M_GT=0.0)
|
||||
|
||||
assert isinstance(result, GVAExtractionResult)
|
||||
assert (
|
||||
abs(result.g_V - true_gV) < 0.15
|
||||
), f"g_V: expected {true_gV}, got {result.g_V}"
|
||||
|
||||
def test_gV_gA_pure_gt(self):
|
||||
"""For pure Gamow-Teller, g_V ≈ 0."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(10, 250, 50)
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
|
||||
# Pure GT: C(W) = g_A² (constant)
|
||||
true_gA = 0.8
|
||||
cw = true_gA**2 + np.random.normal(0, 0.02, len(W))
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
transition_type="F1",
|
||||
use_screening=True,
|
||||
use_exchange=True,
|
||||
use_finite_size=True,
|
||||
use_radiative=True,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(
|
||||
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.02
|
||||
)
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.02,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.02,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
result = extractor.fit_gV_gA(cw_result, M_F=0.0, M_GT=1.0)
|
||||
|
||||
assert isinstance(result, GVAExtractionResult)
|
||||
assert (
|
||||
abs(result.g_A - true_gA) < 0.15
|
||||
), f"g_A: expected {true_gA}, got {result.g_A}"
|
||||
|
||||
def test_gV_gA_mixed_transition(self):
|
||||
"""For mixed transition, extract both g_V and g_A."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(10, 250, 50)
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
|
||||
# Mixed: C(W) = g_V² + g_A² (constant for simplicity)
|
||||
true_gV, true_gA = 0.4, 0.6
|
||||
cw = true_gV**2 + true_gA**2 + np.random.normal(0, 0.02, len(W))
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
transition_type="F1",
|
||||
use_screening=True,
|
||||
use_exchange=True,
|
||||
use_finite_size=True,
|
||||
use_radiative=True,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(
|
||||
config, cw, energies_keV, flux_errors=np.ones_like(cw) * 0.02
|
||||
)
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.02,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.02,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
result = extractor.fit_gV_gA(cw_result, M_F=1.0, M_GT=1.0)
|
||||
|
||||
assert isinstance(result, GVAExtractionResult)
|
||||
# Check that sum of squares is correct
|
||||
reconstructed = result.g_V**2 + result.g_A**2
|
||||
expected = true_gV**2 + true_gA**2
|
||||
assert (
|
||||
abs(reconstructed - expected) < 0.15
|
||||
), f"g_V²+g_A²: expected {expected}, got {reconstructed}"
|
||||
|
||||
def test_gV_gA_summary(self):
|
||||
"""summary() should return a non-empty string with key info."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(10, 250, 50)
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
cw = np.full_like(W, 0.5)
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(config, cw, energies_keV)
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.02,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.02,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
result = extractor.fit_gV_gA(cw_result, M_F=1.0, M_GT=0.0)
|
||||
|
||||
summary = result.summary()
|
||||
assert "g_V" in summary
|
||||
assert "g_A" in summary
|
||||
assert "χ²" in summary
|
||||
assert "p-value" in summary
|
||||
|
||||
|
||||
class TestKurieAnalysis:
|
||||
"""Test Kurie plot analysis."""
|
||||
|
||||
def test_kurie_analysis_allowed_decay(self):
|
||||
"""Kurie analysis on allowed decay should produce reasonable C(W)."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(5, 280, 80)
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
transition_type="F1",
|
||||
use_screening=True,
|
||||
use_exchange=True,
|
||||
use_finite_size=True,
|
||||
use_radiative=True,
|
||||
)
|
||||
spectrum = BetaSpectrum.from_config(config)
|
||||
W, E_MeV = spectrum.get_energy_grid(config)
|
||||
flux_theory = spectrum(W)
|
||||
|
||||
# Interpolate to our energy grid
|
||||
from scipy.interpolate import interp1d
|
||||
|
||||
flux_interp = interp1d(
|
||||
E_MeV * 1000, flux_theory, kind="linear", fill_value=0.0, bounds_error=False
|
||||
)
|
||||
flux = flux_interp(energies_keV)
|
||||
|
||||
# Add noise
|
||||
flux += np.random.normal(0, 0.02 * np.max(flux), len(flux))
|
||||
|
||||
extractor = CWExtractor(config, flux, energies_keV)
|
||||
cw_result, fit_result = extractor.kurie_analysis(endpoint_keV=294.0)
|
||||
|
||||
assert isinstance(cw_result, CWExtractionResult)
|
||||
assert cw_result.fit_result is not None
|
||||
assert np.all(cw_result.cw_values > 0)
|
||||
|
||||
def test_kurie_plot_values(self):
|
||||
"""Kurie plot should produce valid values."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(10, 250, 50)
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
cw = np.ones_like(W) * 1.0
|
||||
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.05,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.05,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
E_plot, kurie = cw_result.kurie_plot()
|
||||
|
||||
assert len(E_plot) == len(energies_keV)
|
||||
assert len(kurie) == len(energies_keV)
|
||||
assert np.all(kurie >= 0), "Kurie values should be non-negative"
|
||||
|
||||
|
||||
class TestEdgeCases:
|
||||
"""Test edge cases and error handling."""
|
||||
|
||||
def test_not_enough_points_for_fit(self):
|
||||
"""Fitting with too few points should raise ValueError."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.array([10.0, 20.0, 30.0])
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
cw = np.ones_like(W)
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(config, cw, energies_keV)
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.05,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.05,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
with pytest.raises(ValueError, match="Not enough data points"):
|
||||
extractor.fit_CW_parametrization(cw_result, parametrization="quadratic")
|
||||
|
||||
def test_not_enough_points_for_gV_gA(self):
|
||||
"""g_V/g_A fit with too few points should raise ValueError."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.array([10.0, 20.0, 30.0])
|
||||
W = np.array([T_to_W(e / 1000.0) for e in energies_keV])
|
||||
cw = np.ones_like(W)
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(config, cw, energies_keV)
|
||||
cw_result = CWExtractionResult(
|
||||
energies_W=W.copy(),
|
||||
energies_keV=energies_keV.copy(),
|
||||
cw_values=cw.copy(),
|
||||
cw_errors=np.ones_like(cw) * 0.05,
|
||||
flux=cw.copy(),
|
||||
flux_errors=np.ones_like(cw) * 0.05,
|
||||
endpoint_W=T_to_W(0.294),
|
||||
endpoint_keV=294.0,
|
||||
)
|
||||
|
||||
with pytest.raises(ValueError, match="Not enough data points"):
|
||||
extractor.fit_gV_gA(cw_result, M_F=1.0, M_GT=1.0)
|
||||
|
||||
def test_default_uncertainties(self):
|
||||
"""Default uncertainties should be sqrt(flux)."""
|
||||
np.random.seed(42)
|
||||
energies_keV = np.linspace(10, 250, 30)
|
||||
flux = np.array(
|
||||
[
|
||||
100.0,
|
||||
200.0,
|
||||
150.0,
|
||||
180.0,
|
||||
120.0,
|
||||
90.0,
|
||||
80.0,
|
||||
70.0,
|
||||
60.0,
|
||||
50.0,
|
||||
45.0,
|
||||
40.0,
|
||||
35.0,
|
||||
30.0,
|
||||
25.0,
|
||||
20.0,
|
||||
18.0,
|
||||
15.0,
|
||||
12.0,
|
||||
10.0,
|
||||
8.0,
|
||||
7.0,
|
||||
6.0,
|
||||
5.0,
|
||||
4.0,
|
||||
3.0,
|
||||
2.0,
|
||||
1.5,
|
||||
1.0,
|
||||
0.5,
|
||||
]
|
||||
)
|
||||
|
||||
config = SpectrumConfig(
|
||||
Z_parent=43,
|
||||
Z_daughter=44,
|
||||
A_number=99,
|
||||
endpoint_MeV=0.294,
|
||||
)
|
||||
|
||||
extractor = CWExtractor(config, flux, energies_keV)
|
||||
expected_unc = np.sqrt(np.maximum(flux, 1.0))
|
||||
np.testing.assert_allclose(extractor.flux_errors, expected_unc)
|
||||
@@ -18,7 +18,7 @@ We test:
|
||||
|
||||
import numpy as np
|
||||
|
||||
from beta_spectrum import CurveFitter, FitConfig
|
||||
from fitter.fit_engine import CurveFitter, FitConfig
|
||||
|
||||
|
||||
class TestBasicFitting:
|
||||
|
||||
@@ -30,10 +30,10 @@ import pytest
|
||||
|
||||
from beta_spectrum import (
|
||||
BetaSpectrum,
|
||||
BetaSpectrumAnalyzer,
|
||||
BranchConfig,
|
||||
SpectrumConfig,
|
||||
)
|
||||
from beta_spectrum.visualize import BetaSpectrumAnalyzer
|
||||
from beta_spectrum.utils import T_to_W
|
||||
|
||||
# ---------------------------------------------------------------------------
|
||||
|
||||
Reference in New Issue
Block a user