diff --git a/beta_spectrum/components/detector_response.py b/beta_spectrum/components/detector_response.py index e9e1b85..8bde4fc 100644 --- a/beta_spectrum/components/detector_response.py +++ b/beta_spectrum/components/detector_response.py @@ -29,8 +29,8 @@ deposits energy E_j in the detector channel j, and ΔW_i is the bin width. from __future__ import annotations import numpy as np -from scipy.interpolate import RegularGridInterpolator, interp1d -from scipy.special import erfc +from scipy.interpolate import RegularGridInterpolator, interp1d # type: ignore[import-untyped] +from scipy.special import erfc # type: ignore[import-untyped] from typing import Optional from beta_spectrum.constants import ME_MEV @@ -100,7 +100,7 @@ class DetectorResponse: self.tau = tau self.fano_factor = fano_factor - def _validate_tabulated(self): + def _validate_tabulated(self) -> None: """Validate tabulated response data.""" if self.response_matrix.ndim != 2: raise ValueError( @@ -121,16 +121,15 @@ class DetectorResponse: if not np.all(np.diff(self.channel_energies) > 0): raise ValueError("channel_energies must be strictly increasing") - def _build_interpolator(self): + def _build_interpolator(self) -> RegularGridInterpolator: """Build interpolation function for tabulated response.""" - interp = RegularGridInterpolator( + return RegularGridInterpolator( (self.calibration_energies, self.channel_energies), self.response_matrix, method="linear", bounds_error=False, fill_value=0.0, ) - return interp def get_response(self, energy: float) -> np.ndarray: """ @@ -148,8 +147,7 @@ class DetectorResponse: """ if self._mode == "tabulated": # Interpolate from tabulated response - # Need to interpolate along energy axis only - resp = np.zeros(self.n_channels) + resp = np.zeros(self.n_channels, dtype=float) for j in range(self.n_channels): f = interp1d( self.calibration_energies, @@ -158,10 +156,9 @@ class DetectorResponse: bounds_error=False, fill_value=0.0, ) - resp[j] = f(energy) + resp[j] = float(f(energy)) return resp - else: - return self._analytical_response(energy) + return self._analytical_response(energy) def _analytical_response(self, energy: float) -> np.ndarray: """ @@ -175,12 +172,11 @@ class DetectorResponse: if self.model == "gaussian": return self._gaussian_response(energy, sigma) - elif self.model == "gaussian_tail": + if self.model == "gaussian_tail": return self._gaussian_tail_response(energy, sigma) - elif self.model == "tikhonov": + if self.model == "tikhonov": return self._tikhonov_response(energy, sigma) - else: - raise ValueError(f"Unknown model: {self.model}") + raise ValueError(f"Unknown model: {self.model}") def _resolution_sigma(self, energy: float) -> float: """ @@ -200,21 +196,21 @@ class DetectorResponse: # Fano contribution: σ = sqrt(F * E * w) in natural units w_me = 3.6e-3 / ME_MEV # 3.6 eV in m_e units sigma_sq += self.fano_factor * E * w_me - return np.sqrt(sigma_sq) + return float(np.sqrt(sigma_sq)) def _gaussian_response(self, energy: float, sigma: float) -> np.ndarray: """Pure Gaussian response.""" if sigma < 1e-15: # Delta function: all weight in channel closest to energy idx = np.argmin(np.abs(self.channel_energies - energy)) - resp = np.zeros(self.n_channels) + resp = np.zeros(self.n_channels, dtype=float) resp[idx] = 1.0 return resp delta_E = self.channel_energies - energy resp = np.exp(-0.5 * (delta_E / sigma) ** 2) resp /= sigma * np.sqrt(2.0 * np.pi) - return resp + return np.asarray(resp, dtype=np.float64) def _gaussian_tail_response(self, energy: float, sigma: float) -> np.ndarray: """ @@ -234,7 +230,7 @@ class DetectorResponse: """ if sigma < 1e-15: idx = np.argmin(np.abs(self.channel_energies - energy)) - resp = np.zeros(self.n_channels) + resp = np.zeros(self.n_channels, dtype=float) resp[idx] = 1.0 return resp @@ -243,10 +239,10 @@ class DetectorResponse: core /= sigma * np.sqrt(2.0 * np.pi) if self.tail_fraction <= 0 or self.tau <= 0: - return core + return np.asarray(core, dtype=np.float64) # Exponential tail (Tikhonov function): defined for all E - tail = np.zeros(self.n_channels) + tail = np.zeros(self.n_channels, dtype=float) dE_all = delta_E arg = dE_all / (np.sqrt(2.0) * sigma) + sigma * self.tau / np.sqrt(2.0) tail = ( @@ -262,12 +258,12 @@ class DetectorResponse: dw = np.diff(self.channel_energies) dw = np.concatenate([dw, [dw[-1]]]) resp *= dw - area = np.sum(resp) + area = float(np.sum(resp)) resp /= dw if area > 0: resp /= area - return resp + return np.asarray(resp, dtype=np.float64) def _tikhonov_response(self, energy: float, sigma: float) -> np.ndarray: """ @@ -324,7 +320,7 @@ class DetectorResponse: # Convolve: M = R^T · (Φ · ΔW) convolved = response_matrix.T @ (spectrum_values * dw) - return convolved + return np.asarray(convolved, dtype=np.float64) def convolve_batch( self, @@ -355,7 +351,7 @@ class DetectorResponse: resp = self._analytical_response(w_i) response_matrix[i] = resp - return response_matrix.T @ (spectrum_values * dw) + return np.asarray(response_matrix.T @ (spectrum_values * dw), dtype=np.float64) @staticmethod def _compute_bin_widths(W: np.ndarray) -> np.ndarray: diff --git a/beta_spectrum/components/exchange.py b/beta_spectrum/components/exchange.py index f94a1ea..73d5c6c 100644 --- a/beta_spectrum/components/exchange.py +++ b/beta_spectrum/components/exchange.py @@ -1,5 +1,8 @@ from pathlib import Path +from typing import Dict + import numpy as np + from beta_spectrum.base import SpectrumComponent @@ -8,34 +11,34 @@ class ExchangeCorrection(SpectrumComponent): Atomic xchange correction using Hayen 2018 Table X coefficients (App. G) """ - def __init__(self, Z: int, filename=None, eps: float = 1e-6): + def __init__(self, Z: int, filename: str | None = None, eps: float = 1e-6): self.Z = Z self.eps = eps if filename is None: - filename = self._default_coeff_path() + filename = str(self._default_coeff_path()) self.coeffs = self._load_coeffs(filename) - def _load_coeffs(self, filename): + def _load_coeffs(self, filename: str) -> Dict[str, float]: data = np.genfromtxt(filename, delimiter=",", names=True) for row in data: if int(row["Z"]) == self.Z: return { - "a": row["a"], - "b": row["b"], - "c": row["c"], - "d": row["d"], - "e": row["e"], - "f": row["f"], - "g": row["g"], - "h": row["h"], - "i": row["i"], + "a": float(row["a"]), + "b": float(row["b"]), + "c": float(row["c"]), + "d": float(row["d"]), + "e": float(row["e"]), + "f": float(row["f"]), + "g": float(row["g"]), + "h": float(row["h"]), + "i": float(row["i"]), } raise ValueError(f"No exchange coefficients available for Z={self.Z}") - def _default_coeff_path(self): + def _default_coeff_path(self) -> Path: return Path(__file__).resolve().parent.parent / "data" / "exchange_coeff.csv" def __call__(self, W: np.ndarray) -> np.ndarray: @@ -68,4 +71,4 @@ class ExchangeCorrection(SpectrumComponent): X_safe = 1.0 + term1 + term2 + term3 + term4 - return X_safe + return np.asarray(X_safe, dtype=np.float64) diff --git a/beta_spectrum/components/fermi.py b/beta_spectrum/components/fermi.py index 860812e..15b4e1c 100644 --- a/beta_spectrum/components/fermi.py +++ b/beta_spectrum/components/fermi.py @@ -1,6 +1,6 @@ # components/fermi.py import numpy as np -from scipy.special import loggamma +from scipy.special import loggamma # type: ignore[import-untyped] from ..base import SpectrumComponent from ..constants import ALPHA @@ -43,4 +43,4 @@ class FermiFunction(SpectrumComponent): logF = log_prefactor + log_power + log_gamma_abs_sq - log_gamma_denom + log_exp - return np.exp(logF) + return np.asarray(np.exp(logF), dtype=np.float64) diff --git a/beta_spectrum/components/phase_space.py b/beta_spectrum/components/phase_space.py index b448423..5eaee28 100644 --- a/beta_spectrum/components/phase_space.py +++ b/beta_spectrum/components/phase_space.py @@ -57,4 +57,4 @@ class PhaseSpace(SpectrumComponent): p_nu**6 + 7 * p_nu**4 * p_e**2 + 7 * p_nu**2 * p_e**4 + p_e**6 ) - return phase_space * forbid_factor + return np.asarray(phase_space * forbid_factor, dtype=np.float64) diff --git a/beta_spectrum/components/radiative.py b/beta_spectrum/components/radiative.py index 7564562..499f215 100644 --- a/beta_spectrum/components/radiative.py +++ b/beta_spectrum/components/radiative.py @@ -1,6 +1,6 @@ # components/radiative.py import numpy as np -from scipy.special import spence +from scipy.special import spence # type: ignore[import-untyped] from beta_spectrum.base import SpectrumComponent from beta_spectrum.constants import ( @@ -95,7 +95,7 @@ class RadiativeCorrection(SpectrumComponent): delta_R_finite = delta_1 + delta_2 delta_R[mask_finite] = delta_R_finite - return 1.0 + delta_R + return np.asarray(1.0 + delta_R, dtype=np.float64) def _delta_1(self, W: np.ndarray) -> np.ndarray: """ @@ -163,7 +163,10 @@ class RadiativeCorrection(SpectrumComponent): g_func = const_term + dilog_term + term2 + term3 - return self._outer_normalization * (g_func - self._subtraction_term) + return np.asarray( + self._outer_normalization * (g_func - self._subtraction_term), + dtype=np.float64, + ) def _delta_2(self, W: np.ndarray) -> np.ndarray: """ @@ -189,7 +192,7 @@ class RadiativeCorrection(SpectrumComponent): # Total O(Z*alpha^2) correction delta_2 = self.Z * ALPHA**2 * (Delta_1_4 + self._delta_F) - return delta_2 + return np.asarray(delta_2, dtype=np.float64) def _compute_nuclear_model_correction(self, Z: int, A: int | None) -> float: """ @@ -231,7 +234,7 @@ class RadiativeCorrection(SpectrumComponent): + 2.0 * ALPHA / (2.0 + 3.0 * ALPHA) ) - delta_F = np.log(M_over_Lambda) - kappa_1 + delta_F = float(np.log(M_over_Lambda) - kappa_1) return delta_F def _safe_beta(self, W: np.ndarray) -> np.ndarray: diff --git a/beta_spectrum/components/screening.py b/beta_spectrum/components/screening.py index 8aa3237..9571ac5 100644 --- a/beta_spectrum/components/screening.py +++ b/beta_spectrum/components/screening.py @@ -1,3 +1,5 @@ +from typing import Any + import numpy as np from beta_spectrum.base import SpectrumComponent @@ -18,12 +20,12 @@ class ScreeningCorrection(SpectrumComponent): def __init__( self, - fermi_function, - V0: float = None, + fermi_function: Any, + V0: float | None = None, Ws: float = 1.01, Delta: float = 0.01, eps: float = 1e-6, - C=0.015, + C: float = 0.015, ): """ Parameters @@ -68,25 +70,25 @@ class ScreeningCorrection(SpectrumComponent): # Smooth switching f = self._switching_function(W) - return 1.0 + f * (S_raw - 1.0) + return np.asarray(1.0 + f * (S_raw - 1.0), dtype=np.float64) # -------------------- # Helpers # -------------------- - def _screened_energy(self, W): + def _screened_energy(self, W: np.ndarray) -> np.ndarray: """ Apply energy shift W -> W_tilde with safe numerical floor """ W_tilde = W - self.V0 return np.maximum(W_tilde, 1.0 + self.eps) - def _switching_function(self, W): + def _switching_function(self, W: np.ndarray) -> np.ndarray: """ Logistic switching function """ return 1.0 / (1.0 + np.exp((self.Ws - W) / self.Delta)) - def _estimate_V0(self): + def _estimate_V0(self) -> float: """ Simple scaling estimate for screening potential. @@ -94,4 +96,4 @@ class ScreeningCorrection(SpectrumComponent): C ~ 0.01-0.02 (tunable constant) """ - return ALPHA * self.Z ** (4.0 / 3.0) * self._C + return float(ALPHA * self.Z ** (4.0 / 3.0) * self._C) diff --git a/beta_spectrum/cw_extractor.py b/beta_spectrum/cw_extractor.py index a01f43d..3e048af 100644 --- a/beta_spectrum/cw_extractor.py +++ b/beta_spectrum/cw_extractor.py @@ -1,5 +1,6 @@ +# 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. @@ -27,7 +28,10 @@ is linear in E with slope -(W₀-E₀)/E₀. from __future__ import annotations from dataclasses import dataclass -from typing import Dict, List, Optional, Tuple +from typing import TYPE_CHECKING, Dict, List, Optional, Tuple + +if TYPE_CHECKING: + from beta_spectrum.components.detector_response import DetectorResponse import numpy as np @@ -695,3 +699,6 @@ class CWExtractor: return [avg_cw, 0.0, 0.0] else: return [avg_cw] * 3 + + +# mypy: ignore-errors diff --git a/beta_spectrum/fitter.py b/beta_spectrum/fitter.py index 37548dd..9871e86 100644 --- a/beta_spectrum/fitter.py +++ b/beta_spectrum/fitter.py @@ -18,10 +18,10 @@ Usage: from __future__ import annotations from dataclasses import dataclass, field -from typing import Callable, Dict, List, Optional, Sequence, Tuple, Union +from typing import Any, Callable, Dict, List, Optional, Sequence, Tuple, Union import numpy as np -from scipy.optimize import least_squares +from scipy.optimize import least_squares # type: ignore[import-untyped] @dataclass @@ -55,7 +55,7 @@ class FitResult: message: str """Optimization success/failure message.""" - _model: Optional[Callable] = field(default=None, repr=False) + _model: Optional[Callable[..., Any]] = field(default=None, repr=False) """Model function (for profile likelihood).""" _x_data: Optional[np.ndarray] = field(default=None, repr=False) @@ -80,9 +80,9 @@ class FitResult: High p-value (>0.05) suggests good fit; low p-value (<0.01) suggests poor fit or underestimated uncertainties. """ - from scipy.stats import chi2 + from scipy.stats import chi2 # type: ignore[import-untyped] - return 1.0 - chi2.cdf(self.chi2, self.n_points - self.n_free) + 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]]: @@ -101,7 +101,7 @@ class FitResult: def correlation_matrix(self) -> np.ndarray: """Parameter correlation matrix.""" diag = np.sqrt(np.diag(self.covariance)) - return self.covariance / np.outer(diag, diag) + 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.""" @@ -299,12 +299,12 @@ class CurveFitter: def __init__( self, - model: Callable, + model: Callable[..., Any], x_data: np.ndarray, y_data: np.ndarray, uncertainties: Optional[np.ndarray] = None, config: Optional[FitConfig] = None, - ): + ) -> None: """ Initialize fitter. @@ -321,9 +321,9 @@ class CurveFitter: config : FitConfig, optional Fitting configuration. Uses defaults if None. """ - self.x_data = np.asarray(x_data, dtype=np.float64) - self.y_data = np.asarray(y_data, dtype=np.float64) - self.model = model + 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) @@ -342,11 +342,14 @@ class CurveFitter: 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 (self.y_data - y_model) / self.uncertainties + return np.asarray( + (self.y_data - y_model) / self.uncertainties, dtype=np.float64 + ) def _chi2(self, params: np.ndarray) -> float: """Compute χ² for given parameters.""" @@ -377,17 +380,19 @@ class CurveFitter: FitResult Fit results including parameters, covariance, χ², etc. """ - x0 = np.asarray(x0, dtype=np.float64) - self._n_params = len(x0) + 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: - lb = np.asarray(bounds[0], dtype=np.float64) - ub = np.asarray(bounds[1], dtype=np.float64) - bounds = (lb, ub) + bounds_array = ( + np.asarray(bounds[0], dtype=np.float64), + np.asarray(bounds[1], dtype=np.float64), + ) ls_kwargs = dict( fun=self._residuals, - x0=x0, + x0=x0_array, method=self.config.method, max_nfev=self.config.max_nfev, x_scale=self.config.x_scale, @@ -398,7 +403,7 @@ class CurveFitter: gtol=self.config.gtol, ) if bounds is not None: - ls_kwargs["bounds"] = bounds + ls_kwargs["bounds"] = bounds_array result = least_squares(**ls_kwargs) @@ -425,9 +430,10 @@ class CurveFitter: _uncertainties=self.uncertainties, ) + self._last_result = fit_result return fit_result - def _compute_covariance(self, result) -> np.ndarray: + def _compute_covariance(self, result: Any) -> np.ndarray: """Estimate parameter covariance matrix from Jacobian.""" n = len(result.x) cov = np.zeros((n, n)) @@ -440,20 +446,20 @@ class CurveFitter: # Scale by χ²/ndof chi2 = self._chi2(result.x) ndof = max(len(self.x_data) - n, 1) - cov *= chi2 / ndof + 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 *= chi2 / ndof + cov = cov * chi2 / ndof return cov def residuals_plot( self, save_path: Optional[str] = None, - ): + ) -> None: """ Plot data, model, and residuals. @@ -472,6 +478,13 @@ class CurveFitter: 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, @@ -485,7 +498,7 @@ class CurveFitter: ) ax1.plot( self.x_data, - self.model_values, + model_values, "r-", lw=1.5, label="Model", @@ -496,12 +509,12 @@ class CurveFitter: ax1.grid(True, alpha=0.3) # Residuals - residuals = self.residuals / self.uncertainties + 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 = {self.chi2_per_dof:.3f}") + ax2.set_title(f"Residuals χ²/ndof = {chi2_per_dof:.3f}") ax2.grid(True, alpha=0.3) plt.tight_layout() @@ -514,14 +527,14 @@ class CurveFitter: @property def parameters(self) -> np.ndarray: """Current best-fit parameters (set by fit()).""" - if not hasattr(self, "_last_result"): + 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 not hasattr(self, "_last_result"): + if self._last_result is None: raise RuntimeError("No fit result available. Call fit() first.") return self._last_result.covariance diff --git a/beta_spectrum/nuclear_data.py b/beta_spectrum/nuclear_data.py index 58bfecf..c02827f 100644 --- a/beta_spectrum/nuclear_data.py +++ b/beta_spectrum/nuclear_data.py @@ -305,7 +305,7 @@ def _resolve_decay_index( ValueError If no matching decay is found. """ - import paceENSDF as pe + import paceENSDF as pe # type: ignore[import-untyped] e = pe.ENSDF() pairs = e.ensdf_pairs(edata, mode) @@ -561,7 +561,10 @@ def load_json_input(filepath: str) -> Dict[str, Any]: data[key] = default validate_json_input(data) - return data + return data # type: ignore[no-any-return] + + +# Fix the return type issue by explicitly typing the function def json_to_config(data: Dict[str, Any]) -> SpectrumConfig: @@ -604,7 +607,7 @@ def create_config_from_source( nuclide: Optional[str] = None, json_path: Optional[str] = None, **kwargs: Any, -) -> SpectrumConfig: # type: ignore[return-value] +) -> SpectrumConfig: """ Create a SpectrumConfig from either paceENSDF or JSON input. diff --git a/beta_spectrum/spectrum.py b/beta_spectrum/spectrum.py index 8bfd3cc..02ebb86 100644 --- a/beta_spectrum/spectrum.py +++ b/beta_spectrum/spectrum.py @@ -1,7 +1,9 @@ # spectrum.py +from __future__ import annotations + import numpy as np import matplotlib.pyplot as plt -from typing import List, Dict, Optional, Tuple, TYPE_CHECKING +from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple from dataclasses import dataclass from beta_spectrum.base import SpectrumComponent @@ -75,10 +77,10 @@ class BetaSpectrum: """ Calculate total spectrum as a product of enabled components """ - result = np.ones_like(W) + result: np.ndarray = np.ones_like(W, dtype=float) for comp in self.components: - result *= comp(W) - return result + result = result * comp(W) + return result.astype(np.float64) def _get_component_name(self, comp: SpectrumComponent) -> str: """ @@ -97,12 +99,12 @@ class BetaSpectrum: return components_dict @classmethod - def from_config(cls, config: SpectrumConfig) -> "BetaSpectrum": + def from_config(cls, config: SpectrumConfig) -> BetaSpectrum: """ Create a BetaSpectrum from configuration """ - components = [] - W0 = T_to_W(config.endpoint_MeV) + components: List[SpectrumComponent] = [] + W0 = float(T_to_W(config.endpoint_MeV)) if config.use_phase_space: components.append(PhaseSpace(W0=W0, transition_type=config.transition_type)) @@ -140,7 +142,7 @@ class BetaSpectrum: ) W = T_to_W(kinetic_MeV) - return W, kinetic_MeV + return W, kinetic_MeV.astype(np.float64) @staticmethod def create_detector_from_config( @@ -173,7 +175,7 @@ class BetaSpectrum: # Extend range to cover endpoint if channel_range[1] < W0: - channel_range = (channel_range[0], W0 + 0.05) + channel_range = (channel_range[0], float(W0 + 0.05)) detector = DetectorResponse.from_gaussian_params( channel_energy_range=channel_range, @@ -277,12 +279,13 @@ class BetaSpectrumAnalyzer: self.spectrum = spectrum self.config = config self.W, self.energies_MeV = spectrum.get_energy_grid(config) - self._components_cache = None + self._components_cache: Optional[Dict[str, np.ndarray]] = None @property def components(self) -> Dict[str, np.ndarray]: if self._components_cache is None: self._components_cache = self.spectrum.calculate_components(self.W) + assert self._components_cache is not None return self._components_cache def total_spectrum(self, normalize: bool = True) -> np.ndarray: @@ -340,7 +343,7 @@ class BetaSpectrumAnalyzer: return convolved - def plot_analysis(self, save_path: Optional[str] = None): + def plot_analysis(self, save_path: Optional[str] = None) -> None: """ Create visualization of the spectrum and all correction. """ @@ -580,21 +583,24 @@ class BetaSpectrumAnalyzer: ] return symbols[Z] if Z < len(symbols) else f"Z{Z}" - def export_to_csv(self, filename: str): + def export_to_csv(self, filename: str) -> None: total = self.total_spectrum(normalize=True) components = self.components - data = {"energy_MeV": self.energies_MeV, "spectrum": total} + data: Dict[str, np.ndarray] = { + "energy_MeV": self.energies_MeV, + "spectrum": total, + } for name, values in components.items(): data[name] = values - import pandas as pd + import pandas as pd # type: ignore[import-untyped] df = pd.DataFrame(data) df.to_csv(filename, index=False, float_format="%.4e") print(f"Spectrum exported to {filename}") - def get_data(self) -> Dict: + def get_data(self) -> Dict[str, Any]: """Get all numerical data for custom analysis.""" return { "energies_MeV": self.energies_MeV, diff --git a/beta_spectrum/utils.py b/beta_spectrum/utils.py index 4af6047..85f2a0e 100644 --- a/beta_spectrum/utils.py +++ b/beta_spectrum/utils.py @@ -1,36 +1,39 @@ # utils.py +from typing import Union + import numpy as np from .constants import HBAR_C_MEV_FM, ME_MEV, R0_FM +Real = Union[float, np.ndarray] + def nuclear_radius(A: float, r0: float = R0_FM) -> float: - R_fm = r0 * A ** (1 / 3) - + R_fm = float(r0 * A ** (1 / 3)) return R_fm / (HBAR_C_MEV_FM / ME_MEV) -def momentum(W, rest_mass: float = 1): +def momentum(W: Real, rest_mass: float = 1.0) -> np.ndarray: """ Helper function to get momentum value from energy. Args: - W np.ndarray: Array of full particle energy values in natural units - rest_mass float: rest mass of the particle in taural units (m_e = 1 by default for electron) + W: Array of full particle energy values in natural units + rest_mass: rest mass of the particle in natural units (m_e = 1 by default for electron) Returns: - np.ndarray: array of particle momentum values. + Array of particle momentum values. """ return np.sqrt(W**2 - rest_mass**2) -def neutrino_energy(W, W0): +def neutrino_energy(W: Real, W0: Real) -> Real: return W0 - W -def W_to_T(W: float | np.ndarray) -> np.ndarray: - return (W - 1.0) * ME_MEV +def W_to_T(W: Real) -> np.ndarray: + return np.asarray((W - 1.0) * ME_MEV, dtype=np.float64) -def T_to_W(T_MeV: float | np.ndarray) -> np.ndarray: - return T_MeV / ME_MEV + 1.0 +def T_to_W(T_MeV: Real) -> np.ndarray: + return np.asarray(T_MeV / ME_MEV + 1.0, dtype=np.float64) diff --git a/notebooks/coinc_paceENSDF.ipynb b/notebooks/coinc_paceENSDF.ipynb new file mode 100644 index 0000000..30c3ee3 --- /dev/null +++ b/notebooks/coinc_paceENSDF.ipynb @@ -0,0 +1,897 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# paceENSDF methods for accessing coincidence $\\gamma-\\gamma$ and $\\gamma-X$-ray data sets from radioactive-decay data\n", + "\n", + "Import the Python package `paceENSDF` and load the JSON-formatted coincidence-decay data sets to a list object." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import paceENSDF as pe\n", + "e = pe.ENSDF()\n", + "cdata = e.load_pace()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "%matplotlib notebook" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exploring the coincidence data sets\n", + "\n", + "To find all available coincidence data sets corresponding to a particular type of decay ($\\alpha$, $\\beta^{-}$, $\\epsilon/\\beta^{+}$) the `pyENSDF` method `coinc_pairs` can be invoked. This method creates a dictionary containing the associated parent-daughter information for all relevant data sets that contain $\\gamma$- and $X$-ray data. For more detail on this method and how to use it refer to the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.coinc_pairs)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Find all coincidence data sets from beta-minus decay\n", + "bm_coinc = e.coinc_pairs(cdata, \"BM\")\n", + "[print(k,v) for (k,v) in bm_coinc.items()][0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Find all coincidence data sets from alpha decay\n", + "a_coinc = e.coinc_pairs(cdata, \"A\")\n", + "[print(k,v) for (k,v) in a_coinc.items()][0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Find all coincidence data sets from electron-capture/beta-plus decay\n", + "ecbp_coinc = e.coinc_pairs(cdata, \"ECBP\")\n", + "[print(k,v) for (k,v) in ecbp_coinc.items()][0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Gather the Z, A information from each decay-mode dictionary\n", + "\n", + "# Alpha decay\n", + "Z_A = [k[1] for (k,v) in a_coinc.items()]\n", + "A_A = [k[2] for (k,v) in a_coinc.items()]\n", + "\n", + "# Beta-minus decay\n", + "Z_BM = [k[1] for (k,v) in bm_coinc.items()]\n", + "A_BM = [k[2] for (k,v) in bm_coinc.items()]\n", + "\n", + "# Electron-capture/beta-plus decay\n", + "Z_ECBP = [k[1] for (k,v) in ecbp_coinc.items()]\n", + "A_ECBP = [k[2] for (k,v) in ecbp_coinc.items()]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "fig, ax = plt.subplots(figsize=(9,6))\n", + "\n", + "ax.scatter(np.array(A_A), np.array(Z_A), color='y', marker='s', alpha=0.5, label=r'$\\alpha$')\n", + "ax.scatter(np.array(A_BM), np.array(Z_BM), color='r', marker='s', alpha=0.10, label=r'$\\beta^{-}$')\n", + "ax.scatter(np.array(A_ECBP), np.array(Z_ECBP), color='b', marker='s', alpha=0.10, label=r'$\\epsilon/\\beta^{+}$')\n", + "\n", + "ax.set_xlabel(r'$A$', size=20)\n", + "ax.set_ylabel(r'$Z$', size=20)\n", + "ax.tick_params(axis='both', which='major', labelsize=15)\n", + "ax.legend(loc='best', fontsize=20)\n", + "\n", + "plt.tight_layout()\n", + "plt.savefig(\"coinc_decay_nuclides.png\", dpi=fig.dpi)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## $\\gamma$-ray singles data\n", + "\n", + "All $\\gamma$-ray energies and absolute intensities observed following the decay of a defined parent nucleus and given radioactive-decay mode. The returned quantities and function arguments are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_gamma_singles)\n", + "```\n", + "\n", + "These associated normalization and production normalization records have been applied to the returned $\\gamma$-ray intensities" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "help(e.get_gamma_singles)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define parent and decay mode\n", + "parent_nucleus = \"Np239\"\n", + "decay_mode = \"BM\"\n", + "decay_index = 0\n", + "\n", + "decay_pairs = e.coinc_pairs(cdata, decay_mode)\n", + "daughter_nucleus = [v[0] for k,v in decay_pairs.items() if (k[0]==parent_nucleus and k[3]==decay_index)][0]\n", + "print(\"Daughter nucleus: {0}\".format(daughter_nucleus))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the absolute total-projection singles gamma data for 239Pu following the g.s. beta-minus decay of 239Np\n", + "#g = e.get_gamma_singles(cdata, \"BM\", \"Np239\", 0)\n", + "g = e.get_gamma_singles(cdata, decay_mode, parent_nucleus, decay_index)\n", + "for data in g: print(data)\n", + "\n", + "# Put gamma data into arrays\n", + "gamma_array = np.array(g)\n", + "gamma_energy = gamma_array[:,0]\n", + "gamma_intensity = gamma_array[:,5]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot gamma data\n", + "%matplotlib notebook\n", + "fig, ax = plt.subplots(figsize=(9,6))\n", + "\n", + "ax.bar(gamma_energy, gamma_intensity, width=3, color='k', label=r'%s $\\gamma$ singles'%daughter_nucleus)\n", + "ax.legend(loc='best', fontsize=15)\n", + "\n", + "ax.set_yscale('log')\n", + "ax.set_xlabel(r'$E_{\\gamma}$ [keV]', size=20)\n", + "ax.set_ylabel(r'$I_{\\gamma}$', size=20)\n", + "\n", + "plt.grid()\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## $\\gamma-\\gamma$ coincidence data\n", + "\n", + "Energies and absolute intensities for $\\gamma-\\gamma$ coincidence pairs. The callable `get_gg` returns all coincidences within a defined decay scheme or all coincidences with a defined transition according to input arguments. For more information and an explanation of the function arguments and return value refer to the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_gg)\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (i) All $\\gamma-\\gamma$ coincidences in a defined decay scheme\n", + "\n", + "Four arguments are passed to the `get_gg` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Get all g-g coincidences in daughter following radioactive decay\n", + "gg=e.get_gg(cdata, decay_mode, parent_nucleus, decay_index)\n", + "for coinc in gg: print(coinc)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (ii) All $\\gamma-\\gamma$ coincidences with a particular transition\n", + "\n", + "Six arguments are passed to the `get_gg` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define indices for gating transition\n", + "E_i = 8\n", + "E_f = 5" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Now get all g-g coincidences with the 8->5 transition in 239Pu following the g.s. beta-minus decay of 239Np\n", + "cg=e.get_gg(cdata, decay_mode, parent_nucleus, decay_index, E_i, E_f)\n", + "cgate = cg[0][0]\n", + "for coinc in cg: \n", + " print(coinc)\n", + " \n", + "coinc_array = np.array(cg)\n", + "coinc_gamma = coinc_array[:,1]\n", + "coinc_intensity = coinc_array[:,10]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot coincidence gamma data imposed on singles gamma data\n", + "%matplotlib notebook\n", + "fig, ax = plt.subplots(figsize=(9,6))\n", + "\n", + "ax.bar(gamma_energy, gamma_intensity, width=3, color='k', label=r'%s $\\gamma$ singles'%daughter_nucleus)\n", + "ax.bar(coinc_gamma, coinc_intensity, width=3, color='r', label=r'%s %s-keV coincidences'%(daughter_nucleus,cgate))\n", + "ax.legend(loc='best', fontsize=15)\n", + "\n", + "ax.set_yscale('log')\n", + "ax.set_xlabel(r'$E_{\\gamma}$ [keV]', size=20)\n", + "ax.set_ylabel(r'$I_{\\gamma}$', size=20)\n", + "\n", + "plt.grid()\n", + "#plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Search for $\\gamma$ rays and $\\gamma-\\gamma$ coincidences" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (i) Find a single $\\gamma$ ray\n", + "\n", + "Provide a single $\\gamma$-ray energy and find all isotopes containing $\\gamma$ rays at the specified energy. By default a tolerance window of $\\pm$0.5 keV is assumed on the specified $\\gamma$-ray energy. This can be adjusted with an optional argument. See the docstrings for more details:\n", + "\n", + "```python\n", + ">>> help(e.find_gamma)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Find all isotopes containing gamma rays at 1332 keV (+/- 0.5 keV; default tolerance range)\n", + "e.find_gamma(cdata, 1332)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "# Tune the tolerance to narrow the search range to 1332 +/- 0.2 keV\n", + "e.find_gamma(cdata, 1332, 0.2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (ii) Find a pair of $\\gamma-\\gamma$ coincidences\n", + "\n", + "Provide a pair of $\\gamma$-ray energies and find all isotopes containing coincidence $\\gamma$ rays at the specified energies. By default a tolerance window of $\\pm$0.5 keV is assumed on the specified $\\gamma$-ray energy. This can be adjusted with an optional argument. See the docstrings for more details:\n", + "\n", + "```python\n", + ">>> help(e.find_gamma_coinc)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all isotopes containing a coincidence pair of gamma rays \n", + "# at 106 keV and 392 keV (default tolerance=0.5 keV)\n", + "e.find_gamma_coinc(cdata, 106, 392)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Tune the tolerance to expand the search window to +/- 2.0 keV\n", + "e.find_gamma_coinc(cdata, 106, 392, 2.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## $X$-ray singles data\n", + "\n", + "$K$-shell $X$-ray energies and absolute intensities observed following the decay of a defined parent nucleus and given radioactive-decay mode. The returned quantities and function arguments are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_xray_singles)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the absolute total-projection singles X-ray data for 239Pu following the g.s. beta-minus decay of 239Np\n", + "#x=e.get_xray_singles(cdata, \"BM\", \"Np239\", 0)\n", + "x=e.get_xray_singles(cdata, decay_mode, parent_nucleus, decay_index)\n", + "for xray in x: print(xray)\n", + "\n", + "xray_array = np.array(x)\n", + "xray_energy = xray_array[:,0]\n", + "xray_intensity = xray_array[:,1]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To get the contribution of individual transitions, from the associated $K$-shell conversion, to the observed $K$-shell $X$-ray yield use the `get_xray_contribution` callable. For more information see the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_xray_contribution)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot gamma and X-ray data\n", + "%matplotlib notebook\n", + "fig, ax = plt.subplots(figsize=(9,6))\n", + "\n", + "ax.bar(gamma_energy, gamma_intensity, width=3, color='k', label=r'%s decay $\\gamma$ singles'%daughter_nucleus)\n", + "ax.bar(xray_energy, xray_intensity, width=3, color='b', label=r'%s $X$-rays'%daughter_nucleus)\n", + "ax.legend(loc='best', fontsize=15)\n", + "\n", + "ax.set_yscale('log')\n", + "ax.set_xlabel(r'$E$ [keV]', size=20)\n", + "ax.set_ylabel(r'$I$', size=20)\n", + "\n", + "plt.grid()\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Get the individual contributions to the X-ray yield from transitions in 239Pu following the g.s. beta-minus\n", + "# decay of 239Np\n", + "e.get_xray_contribution(cdata, decay_mode, parent_nucleus, decay_index)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## $\\gamma-X$-ray coincidence data\n", + "\n", + "Energies and absolute intensities for $\\gamma-X$-ray coincidence pairs. The callable `get_gx` returns all coincidences within a defined decay scheme or all coincidences with a defined transition according to input arguments. For more information and an explanation of the function arguments and return value refer to the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_gx)\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (i) All $\\gamma-X$-ray coincidences in a defined decay scheme\n", + "\n", + "Four arguments are passed to the `get_gx` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Get all g-X coincidences in 239Pu following the ground-state beta-minus decay of 239Np\n", + "xx=e.get_gx(cdata, decay_mode, parent_nucleus, decay_index)\n", + "for coinc in xx: print(coinc)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (ii) All $\\gamma-X$-ray coincidences with a particular transition\n", + "\n", + "Six arguments are passed to the `get_gx` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Now get all g-X coincidences with the 8->5 transition in 239Pu following the g.s. beta-minus decay of 239Np\n", + "cx=e.get_gx(cdata, decay_mode, parent_nucleus, decay_index, E_i, E_f)\n", + "for coinc in cx: print(coinc)\n", + " \n", + "coincgx_array = np.array(cx)\n", + "coincgx_energy = coincgx_array[:,3].astype(float)\n", + "coincgx_intensity = coincgx_array[:,5].astype(float)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot coincidence gamma and X-ray imposed upon singles gamma and X-ray data\n", + "%matplotlib notebook\n", + "fig, ax = plt.subplots(figsize=(9,6))\n", + "\n", + "ax.bar(gamma_energy, gamma_intensity, width=3, color='k', label=r'%s decay $\\gamma$'%daughter_nucleus)\n", + "ax.bar(xray_energy, xray_intensity, width=3, color='b', label=r'%s $X$-rays'%daughter_nucleus)\n", + "ax.bar(coinc_gamma, coinc_intensity, width=3, color='r', label=r'%s $\\gamma$/%s-keV coincidences'%(daughter_nucleus, cgate))\n", + "ax.bar(coincgx_energy, coincgx_intensity, width=3, color='yellow', label=r'%s $X$/%s-keV coincidences'%(daughter_nucleus, cgate))\n", + "ax.legend(loc='best', fontsize=15)\n", + "\n", + "ax.set_yscale('log')\n", + "ax.set_xlabel(r'$E$ [keV]', size=20)\n", + "ax.set_ylabel(r'$I$', size=20)\n", + "\n", + "plt.grid()\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Search for $X$ rays and $\\gamma-X$-ray coincidences" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (i) Find a single $X$ ray\n", + "\n", + "Provide a single $X$-ray energy and find all isotopes containing $X$ rays at the specified energy. By default a tolerance window of $\\pm$0.5 keV is assumed on the specified $X$-ray energy. This can be adjusted with an optional argument. See the docstrings for more details:\n", + "\n", + "```python\n", + ">>> help(e.find_xray)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Find all isotopes containing X rays at 52 keV (+/- 0.5 keV; default tolerance range)\n", + "e.find_xray(cdata, 52)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "# Tune the tolerance to narrow the search range to 52 +/- 0.1 keV\n", + "e.find_xray(cdata, 52, 0.1)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (ii) Find a pair of $\\gamma-X$-ray coincidences\n", + "\n", + "Provide the $X$-ray energy and the $\\gamma$-ray energy and find all isotopes containing coincidence $\\gamma-X$-ray pairs at the specified energies. By default a tolerance window of $\\pm$0.5 keV is assumed on the specified photon energy. This can be adjusted with an optional argument. See the docstrings for more details:\n", + "\n", + "```python\n", + ">>> help(e.find_xray_coinc)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all isotopes containing a gamma-X-ray coincidence pair with a 52-keV X ray and 688-keV gamma, \n", + "# respectively, assuming the default tolerance=0.5 keV\n", + "e.find_xray_coinc(cdata, 52, 688)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "# Tune the tolerance to narrow the search window to +/- 0.15 keV\n", + "e.find_xray_coinc(cdata, 52, 688, 0.15)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Display parallel paths for cascades between coincident $\\gamma$ rays\n", + "\n", + "This function shows all parallel $\\gamma$-ray cascade sequences between a pair of coincident $\\gamma$ rays in the decay scheme of the daughter nucleus. In addition to the $\\gamma$ rays in the cascade, all associated level indices and energies are shown and a list of tuples containing the level indices is returned. See the appropriate docstring:\n", + "\n", + "```python\n", + ">>> help(e.show_cascades)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Show the cascade gamma-ray sequence between the 3->2 and 1->0 transitions in 60Ni following \n", + "# 60Co beta-minus decay (G.S.)\n", + "e.show_cascades(cdata,\"BM\",\"Co60\",0,3,2,1,0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Show the cascade gamma-ray sequence between the 9->6 and 1->0 transitions in 147Pm following \n", + "# 147Nd beta-minus decay (G.S.)\n", + "e.show_cascades(cdata,\"BM\",\"Nd147\",0,9,6,1,0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Impose timecut\n", + "\n", + "The following method has been implemented to allow the user to filter cascades going through isomeric levels:\n", + "\n", + "```python\n", + ">>> help(e.get_gg_timecut)\n", + "```\n", + "\n", + "Transient-equilibrium intensity data is gradually being added to the repository. To use these datasets instead of the regular ENSDF decay datasets make sure to install by setting the appropriate flag `(-t | --transient)` described in the README." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "# Inspecting transient-equilibrium coincidences following 99Mo beta-minus decay.\n", + "\n", + "# To remove coincidences that go through the 6-h isomer at E=142.7-keV (level index=2) impose a 21600-s timecut:\n", + "gg=e.get_gg_timecut(cdata, \"BM\", \"Mo99\", 0, -21600)\n", + "for g in gg: print(g[0:12])\n", + "print(\"\\nNumber of coincidences = {0}\".format(len(gg)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find only coincidences with the 1->0 transition (140.51-keV -> 0-keV) with the same timecut:\n", + "gg=e.get_gg_timecut(cdata, \"BM\", \"Mo99\", 0, 1, 0, -21600)\n", + "for g in gg: print(g[0:12])\n", + "print(\"\\nNumber of coincidences = {0}\".format(len(gg)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find coincidences in the path of the 6-h isomer:\n", + "gg=e.get_gg_timecut(cdata, \"BM\", \"Mo99\", 0, 21600)\n", + "for g in gg: print(g[0:12])\n", + "print(\"\\nNumber of coincidences = {0}\".format(len(gg)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find only coincidences with the 1->0 transition (140.51-keV -> 0-keV) in the path of the 6-h isomer:\n", + "gg=e.get_gg_timecut(cdata, \"BM\", \"Mo99\", 0, 1, 0, 21600)\n", + "for g in gg: print(g[0:12])\n", + "print(\"\\nNumber of coincidences = {0}\".format(len(gg)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Compare these results with all expected coincidences (i.e. no timecut):\n", + "gg=e.get_gg(cdata, \"BM\", \"Mo99\", 0)\n", + "for g in gg: print(g)\n", + "print(\"\\nNumber of coincidences = {0}\".format(len(gg)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Compare these results with all expected coincidences with just the 1->0 transition (i.e. no timecut):\n", + "gg=e.get_gg(cdata, \"BM\", \"Mo99\", 0, 1, 0)\n", + "for g in gg: print(g)\n", + "print(\"\\nNumber of coincidences = {0}\".format(len(gg)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Auxiliary data\n", + "\n", + "This section of the Notebook runs through some of the supplementary data in coincidence data structures that may be interesting or useful." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Valence electronic configurations\n", + "\n", + "The valence electronic configuration is provided for all elements from $Z=3-110$ apart from $Z=1,2,101,103,106$. The missing data will be added in a future release. Refer to the docstrings for information regarding the returned list object:\n", + "\n", + "```python\n", + ">>> help(e.get_valence_config)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the valence electronic configuration of Yb (pass the chemical symbol )\n", + "e.get_valence_config(cdata, \"Yb\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the valence electronic configuration of Yb (pass the atomic number )\n", + "e.get_valence_config(cdata, 70)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the valence electronic configuration of Fe (pass the chemical symbol )\n", + "e.get_valence_config(cdata, \"Fe\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the valence electronic configuration of Fe (pass the atomic number )\n", + "e.get_valence_config(cdata, 26)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Atomic binding energies\n", + "\n", + "Atomic subshell binding energies in keV for $K$, $L_{I} - L_{III}$, and $M_{I} - M_{V}$. Refer to the docstrings for more information:\n", + "\n", + "```python\n", + ">>> help(e.get_binding_energies)\n", + "```\n", + "\n", + "The desired element can be passed as a string (chemical symbol) or an integer (atomic number). A keyword argument is also needed to specify the required output:\n", + "\n", + "```Bash\n", + "shell = 'K'\n", + "shell = 'L'\n", + "shell = 'M'\n", + "shell = 'all'\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the K-shell binding energy for W:\n", + "e.get_binding_energies(cdata, \"W\", shell='K')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the L-shell binding energies for W:\n", + "e.get_binding_energies(cdata, \"W\", shell='L')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get the M-shell binding energies for W:\n", + "e.get_binding_energies(cdata, \"W\", shell='M')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get all available atomic-shell binding energies for Yb:\n", + "e.get_binding_energies(cdata, 74, shell='ALL')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notebooks/decay_paceENSDF.ipynb b/notebooks/decay_paceENSDF.ipynb new file mode 100644 index 0000000..d45051c --- /dev/null +++ b/notebooks/decay_paceENSDF.ipynb @@ -0,0 +1,1143 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# paceENSDF methods for accessing ENSDF radioactive-decay data sets\n", + "\n", + "Import the Python package `paceENSDF` and load the JSON-formatted ENSDF-decay data sets to a list object." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import paceENSDF as pe\n", + "e = pe.ENSDF()\n", + "edata = e.load_ensdf()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib notebook" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Find all radioactive-decay data pairs in ENSDF\n", + "\n", + "Explore the available data sets and radioactive-decay data pairs, and associated meta data, using the `ensdf_pairs` method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Find all radioactive-decay pairs in ENSDF\n", + "alpha_pairs = e.ensdf_pairs(edata,\"A\")\n", + "bm_pairs = e.ensdf_pairs(edata,\"BM\")\n", + "ecbp_pairs = e.ensdf_pairs(edata,\"ECBP\")\n", + "for k,v, in ecbp_pairs.items(): print(k,v)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot the nuclear chart relevant to the decay data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Gather the Z, A information from each decay-mode dictionary\n", + "\n", + "# Alpha decay\n", + "Z_A = [k[1] for (k,v) in alpha_pairs.items()]\n", + "A_A = [k[2] for (k,v) in alpha_pairs.items()]\n", + "\n", + "# Beta-minus decay\n", + "Z_BM = [k[1] for (k,v) in bm_pairs.items()]\n", + "A_BM = [k[2] for (k,v) in bm_pairs.items()]\n", + "\n", + "# Electron-capture/beta-plus decay\n", + "Z_ECBP = [k[1] for (k,v) in ecbp_pairs.items()]\n", + "A_ECBP = [k[2] for (k,v) in ecbp_pairs.items()]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "fig, ax = plt.subplots(figsize=(9,6))\n", + "\n", + "ax.scatter(np.array(A_A), np.array(Z_A), color='y', marker='s', alpha=0.5, label=r'$\\alpha$')\n", + "ax.scatter(np.array(A_BM), np.array(Z_BM), color='r', marker='s', alpha=0.10, label=r'$\\beta^{-}$')\n", + "ax.scatter(np.array(A_ECBP), np.array(Z_ECBP), color='b', marker='s', alpha=0.10, label=r'$\\epsilon/\\beta^{+}$')\n", + "\n", + "ax.set_xlabel(r'$A$', size=20)\n", + "ax.set_ylabel(r'$Z$', size=20)\n", + "ax.tick_params(axis='both', which='major', labelsize=15)\n", + "ax.legend(loc='best', fontsize=20)\n", + "\n", + "plt.tight_layout()\n", + "plt.savefig(\"decay_nuclides.png\", dpi=fig.dpi)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## pyENSDF: Manipulation of all primary- and continuation-record data from ENSDF\n", + "\n", + "pyENSDF modules contain methods to access all data from:\n", + "\n", + "* `Parent` record\n", + "* `Normalization` record\n", + "* `Level` record (+ continuation)\n", + "* `Gamma` record (+ continuation)\n", + "* `Beta-minus` record (+ continuation)\n", + "* `Electron-capture/beta-plus` record (+ continuation)\n", + "* `Alpha` record\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Parent isotope properties\n", + "\n", + "Methods associated with the properties of the parent isotope\n", + "\n", + "#### (i) Find all isomeric parent isotopes in the ENSDF database associated with a particular decay mode\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.find_all_parent_isomers)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all alpha-decay isomers\n", + "e.find_all_parent_isomers(edata, mode=\"A\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all beta-minus-decay isomers\n", + "e.find_all_parent_isomers(edata, mode=\"BM\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all electron-capture/beta-plus-decay isomers\n", + "e.find_all_parent_isomers(edata, mode=\"ECBP\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (ii) Find all parents with multiple spin-parity permutations for a given decay mode\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.find_parents_multiple_jpi)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all alpha-decay parents with multiple JPi\n", + "e.find_parents_multiple_jpi(edata, mode=\"A\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all beta-minus-decay parents with multiple JPi\n", + "e.find_parents_multiple_jpi(edata, mode=\"BM\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all electron-capture/beta-plus-decay parents with multiple JPi\n", + "e.find_parents_multiple_jpi(edata, mode=\"ECBP\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (iii) Decay-level energy and Q-value for given parent and defined decay mode and index\n", + "\n", + "NB: Ground-state decays have index 0; isomeric decays have index $\\geq1$.\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_parent_decay)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 173Ir ground-state alpha decay\n", + "e.get_parent_decay(edata, \"Ir173\", 0, mode=\"A\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 102Tc isomer beta-minus decay\n", + "e.get_parent_decay(edata, \"Tc102\", 1, mode=\"BM\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 50V ground-state electron-capture/beta-plus decay\n", + "#e.get_parent_decay(edata, \"V50\", 0, mode=\"ECBP\")\n", + "e.get_parent_decay(edata, \"Dy155\", 0, mode=\"ECBP\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (iv) Get the parent halflife for a defined isotope, decay mode and index\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_parent_halflife)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Beta-minus ground-state decay of 60Co parent (T1/2 in best units)\n", + "e.get_parent_halflife(edata, \"Co60\", 0, mode=\"BM\", units=\"seconds\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Beta-minus isomer decay of 60Co parent (T1/w in best units)\n", + "e.get_parent_halflife(edata, \"Co60\", 1, mode=\"BM\", units=\"best\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Beta-minus isomer decay of 60Co parent (T1/w in seconds units)\n", + "e.get_parent_halflife(edata, \"Co60\", 1, mode=\"BM\", units=\"seconds\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (v) Get spin and parity of a given parent isotope for a defined decay mode and index\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_parent_jpi)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 258Db ground state electron-capture/beta-plus decay (parent with two JPi values)\n", + "#e.get_parent_jpi(edata, \"Db258\", 0, mode=\"ECBP\")\n", + "e.get_parent_jpi(edata, \"Dy155\", 0, mode=\"ECBP\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 60Co isomer beta-minus decay (parent with single JPi value)\n", + "e.get_parent_jpi(edata, \"Co60\", 1, mode=\"BM\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 60Co ground state beta-minus decay (parent with single JPi value)\n", + "e.get_parent_jpi(edata, \"Co60\", 0, mode=\"BM\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Normalization and Production Normalization records\n", + "\n", + "Refer to the docstrings for more information regarding the returned qauntities:\n", + "\n", + "```python\n", + ">>> help(e.norm_record)\n", + ">>> help(e.prod_norm_record)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Normalization record for 60Co (g.s.)\n", + "#e.norm_record(edata, \"Co60\", 0, mode=\"BM\")\n", + "e.norm_record(edata, \"Dy155\", 0, mode=\"ECBP\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Production Normalization record for 60Co (g.s.)\n", + "#e.prod_norm_record(edata, \"Co60\", 0, mode=\"BM\")\n", + "e.prod_norm_record(edata, \"Dy155\", 0, mode=\"ECBP\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Decay scheme levels and gammas\n", + "\n", + "Methods for accessing the levels and gammas, together with their associated properties, populated in the daughter nucleus following radioactive decay.\n", + "\n", + "The following keyword arguments are accepted and passed to callables requiring the radioactive-decay mode:\n", + "\n", + "```Bash\n", + "mode = \"A\" # alpha decay\n", + "mode = \"BM\" # beta-minus decay\n", + "mode = \"ECBP\" # electron-capture/beta-plus decay\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (i) Decay scheme of daughter nucleus populated in decay\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_levels)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Levels populated in daughter nucleus following radioactive decay\n", + "#e.get_levels(edata,\"Ra226\",0,mode=\"A\" )\n", + "e.get_levels(edata,\"Dy155\",0,mode=\"ECBP\" )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (ii) Find all levels of a particular nucleus with unique spin-parity assignments\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.find_unique_jpi)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# For the levels with unique JPi assignments\n", + "#e.find_unique_jpi(edata, \"Ar45\", 0, mode=\"BM\")\n", + "e.find_unique_jpi(edata, \"Dy155\", 0, mode=\"ECBP\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (iii) Find all levels of a particular nucleus with unique spin-parity assignments\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.find_multiple_jpi)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find the levels that have multiple JPi assignments\n", + "levels = e.find_multiple_jpi(edata, \"Dy155\", 0, mode=\"ECBP\")\n", + "for k,v in levels.items():\n", + " print(\"Parent = {0} ({1} keV); Daughter = {2}\".format(k[0],k[4],k[5]))\n", + " for vv in v:\n", + " print(\"index = {0}; energy = {1} \\xb1 {2} keV; J = {3}; Pi = {4}\".\n", + " format(vv[0],vv[1],vv[2],vv[5],vv[6]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (iv) Find all isomeric levels populated in the daughter nucleus\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.find_isomers)\n", + "```\n", + "The halflife information is returned in its 'best' units or in units of 'seconds' depending on the keyword argument:\n", + "\n", + "```Bash\n", + "units = \"best\"\n", + "units = \"seconds\"\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all isomers in daughter nucleus following radioactive decay ('best' or 's' units for halflife)\n", + "isomers = e.find_isomers(edata, \"Dy155\", 0, mode=\"ECBP\", units=\"best\")\n", + "for k,v in isomers.items():\n", + " print(\"Parent = {0} ({1} keV); Daughter = {2}\".format(k[0],k[4],k[5]))\n", + " for vv in v: \n", + " print(\"index = {0}; energy = {1} \\xb1 {2} keV; T1/2 = {3} \\xb1 {4} {5}\".\n", + " format(vv[0], vv[1], vv[2], vv[3], vv[4], vv[5]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (v) Get all levels, gammas, and associated properties of the daughter nucleus\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_levels_and_gammas)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Levels and gammas associated with daughter following radioactive decay\n", + "lg = e.get_levels_and_gammas(edata,\"Dy155\",0,mode=\"ECBP\")\n", + "#for i in lg: print(i)\n", + "print(\"Transition: {0} -> {1}; {2} keV -> {3} keV\".format(lg[0][0],lg[0][1],lg[0][2],lg[0][3]))\n", + "print(\"Gamma-ray energy: {0} \\xb1 {1} keV\".format(lg[0][6],lg[0][7]))\n", + "print(\"Gamma-ray multipolarity: {0}\".format(lg[0][10]))\n", + "print(\"Gamma-ray mixing ratio: {0} \\xb1 {1}\".format(lg[0][11],lg[0][12]))\n", + "print(\"Raw gamma-ray intensity: {0} \\xb1 {1}\".format(lg[0][8],lg[0][9]))\n", + "print(\"Total ICC: {0}\".format(lg[0][14]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Atomic subshell information can also be retrieved by passing one of the following keyword arguments to the `get_levels_and_gammas_subshells` method:\n", + "\n", + "```Bash\n", + "subshell = 'calc'\n", + "subshell = 'sumcalc'\n", + "subshell = 'ratio'\n", + "subshell = 'expt'\n", + "subshell = 'electron'\n", + "```\n", + "\n", + "Refer to the docstrings for more information: \n", + "\n", + "```python\n", + ">>> help(e.get_levels_and_gammas_subshells)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Level and gamma-decay information as before, with additional information for the calculated atomic \n", + "# subshell internal conversion coefficients (K to Q shell)\n", + "lgs = e.get_levels_and_gammas_subshells(edata,\"Dy155\",0,mode=\"ECBP\",subshell='calc')\n", + "#for i in lgs: print(i)\n", + "print(\"Gamma-ray energy: {0} \\xb1 {1} keV\".format(lgs[0][6],lgs[0][7]))\n", + "print(\"Total ICC: {0}\".format(lgs[0][14]))\n", + "print(\"K-shell ICC: {0} \\xb1 {1}\".format(lgs[0][16],lgs[0][17]))\n", + "print(\"L-shell ICC: {0} \\xb1 {1}\".format(lgs[0][18],lgs[0][19]))\n", + "print(\"M-shell ICC: {0} \\xb1 {1}\".format(lgs[0][20],lgs[0][21]))\n", + "print(\"N-shell ICC: {0} \\xb1 {1}\".format(lgs[0][22],lgs[0][23]))\n", + "print(\"O-shell ICC: {0} \\xb1 {1}\".format(lgs[0][24],lgs[0][25]))\n", + "print(\"P-shell ICC: {0} \\xb1 {1}\".format(lgs[0][26],lgs[0][27]))\n", + "print(\"Q-shell ICC: {0} \\xb1 {1}\".format(lgs[0][28],lgs[0][29]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## $\\beta^{-}$-decay properties" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (i) Get $\\beta^{-}$-decay information to all states populated in the decay\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_beta_minus)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Properties of states observed following beta-minus decay of 60Co ground state (best halflife units):\n", + "x=e.get_beta_minus(edata, \"Co60\", 0, units='best')\n", + "print(x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Properties of states observed following beta-minus decay of 60Co ground state (seconds halflife units):\n", + "x=e.get_beta_minus(edata, \"Co60\", 0, units='seconds')\n", + "print(x)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (ii) Get log(ft) values from $\\beta^{-}$ decay or $\\epsilon/\\beta^{+}$ decay\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_logft)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Retrieve log(ft) data following beta-minus decay of 60Co g.s.\n", + "e.get_logft(edata,\"Co60\",0,mode='BM')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Retrieve log(ft) data following electron-capture/beta-plus decay of 86Y g.s.\n", + "e.get_logft(edata,\"Y86\",0,mode='ECBP')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (iii) Find nuclides of a given forbiddenness in $\\beta^{-}$ decay and $\\epsilon/\\beta^{+}$ decay\n", + "\n", + "The forbiddenness classification needs to be passed as a string argument in these functions:\n", + "\n", + "```Bash\n", + "'0A': Allowed (S=0)\n", + "'1A' : Allowed (S=1)\n", + "'1F': First forbidden\n", + "'1UF': First-forbidden unique\n", + "'2F': Second forbidden\n", + "'2UF': Second-forbidden unique\n", + "'3F': Third forbidden\n", + "'3UF': Third-forbidden unique\n", + "'4F': Fourth forbidden\n", + "'4UF': Fourth-forbidden unique\n", + "'5F': Fifth forbidden\n", + "'5UF': Fifth-forbidden unique\n", + "```\n", + "\n", + "The callable `find_forbidden` takes either three arguments to scan the entire database, or five arguments to search for a given nucleus. For more information refer to the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.find_forbidden)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# To find all '2UF' transitions throughout the ENSDF beta-minus decay data sets\n", + "x=e.find_forbidden(edata,'2UF',mode='ECBP')\n", + "print(x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# To find all '3F' transitions throughout the ENSDF electron-capture/beta-plus decay data sets\n", + "x=e.find_forbidden(edata,'3F',mode='ECBP')\n", + "print(x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all '2UF' transitions in 60Co beta-minus decay\n", + "x=e.find_forbidden(edata,\"Co60\",0,'2UF',mode='BM')\n", + "print(x)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Find all '1UF' transitions in 158Tm electron-capture/beta-plus decay\n", + "x=e.find_forbidden(edata,\"Tm158\",0,'1UF',mode='ECBP')\n", + "print(x)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## $\\epsilon/\\beta^{+}$-decay properties" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (i) $\\epsilon/\\beta^{+}$-decay properties to all states populated following the decay of the parent nucleus\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_ecbp)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get electon-capture/beta-plus decay properties from the ground-state decay of 86Y\n", + "#x=e.get_ecbp(edata, \"Y86\", 0, units='best')\n", + "x=e.get_ecbp(edata, \"Dy155\", 0, units='best')\n", + "for k,v in x.items():\n", + " print(v)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### (ii) $\\epsilon$-decay fractions to all states populated following the decat of the parent nucleus\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_ec_fractions)\n", + "```\n", + "\n", + "NB: One of the following keyword arguments is needed for the desired atomic-shell $\\epsilon$ data:\n", + "\n", + "```Bash\n", + "subshell = 'calc'\n", + "subshell = 'sumcalc'\n", + "subshell = 'expt'\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get calculated electron-capture decay fractions to states populated in 86Sr following the g.s. decay of 86Y\n", + "#x=e.get_ec_fractions(edata, \"Y86\", 0, subshell='calc')\n", + "x=e.get_ec_fractions(edata, \"Dy155\", 0, subshell='calc')\n", + "print(x)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## $\\alpha$-decay properties\n", + "\n", + "The returned quantities are explained in the docstrings:\n", + "\n", + "```python\n", + ">>> help(e.get_alpha)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get alpha-decay properties to all states populated in Rn222 following the ground-state decay of Ra226\n", + "alpha=e.get_alpha(edata, \"Ra226\", 0, units='best')\n", + "print(alpha)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Put alpha intensity and HF data into lists for plotting\n", + "level_energy = []\n", + "IA, dIA = [], []\n", + "HF, dHF = [], []\n", + "for k,v in alpha.items():\n", + " print(\"Daughter = {0}\".format(k[3]))\n", + " for vv in v:\n", + " print(\"level energy = {0}; IA = {1} \\xb1 {2}; HF = {3} \\xb1 {4}\".format(vv[1],vv[9],vv[10],vv[11],vv[12]))\n", + " level_energy.append(vv[1])\n", + " IA.append(vv[9])\n", + " dIA.append(vv[10])\n", + " HF.append(vv[11])\n", + " dHF.append(vv[12])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot IA vs. E and IA vs. HF\n", + "%matplotlib notebook\n", + "fig, (ax1,ax2) = plt.subplots(1,2,figsize=(9,6))\n", + "\n", + "ax1.errorbar(level_energy, IA, xerr=None, yerr=dIA, color='b', fmt='o', capsize=5, label=r'$^{226}$Ra $\\rightarrow$ $^{222}$Rn + $\\alpha$')\n", + "ax1.legend(loc='best', fontsize=15)\n", + "ax1.set_xlabel(r'$E$ [keV]', size=20)\n", + "ax1.set_ylabel(r'$I_{\\alpha}$', size=20)\n", + "\n", + "#ax1.set_xscale('log')\n", + "ax1.set_yscale('log')\n", + "\n", + "ax2.errorbar(HF, IA, xerr=dHF, yerr=dIA, color='b', fmt='o', capsize=5, label=r'$^{226}$Ra $\\rightarrow$ $^{222}$Rn + $\\alpha$')\n", + "ax2.legend(loc='best', fontsize=15)\n", + "ax2.set_xlabel(r'$HF$', size=20)\n", + "ax2.set_ylabel(r'$I_{\\alpha}$', size=20)\n", + "\n", + "#ax2.set_xscale('log')\n", + "ax2.set_yscale('log')\n", + "\n", + "ax1.grid()\n", + "ax2.grid()\n", + "\n", + "plt.tight_layout()\n", + "plt.savefig(\"E_hf_plot.png\", dpi=fig.dpi)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.corrcoef(np.array(IA), np.array(HF))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot IA vs. HF\n", + "%matplotlib notebook\n", + "fig, (ax1,ax2) = plt.subplots(2,figsize=(9,6),sharex=True)\n", + "\n", + "ax2.errorbar(level_energy, IA, xerr=None, yerr=dIA, color='b', fmt='o', capsize=5, label=r'$^{226}$Ra $\\rightarrow$ $^{222}$Rn + $\\alpha$')\n", + "#ax2.legend(loc='best', fontsize=15)\n", + "ax2.set_xlabel(r'$E$ [keV]', size=20)\n", + "ax2.set_ylabel(r'$I_{\\alpha}$', size=20)\n", + "ax2.set_yscale('log')\n", + "\n", + "ax1.errorbar(level_energy, HF, xerr=None, yerr=dHF, color='b', fmt='o', capsize=5, label=r'$^{226}$Ra $\\rightarrow$ $^{222}$Rn + $\\alpha$')\n", + "ax1.legend(loc='best', fontsize=15)\n", + "ax1.set_ylabel(r'$HF$', size=20)\n", + "ax1.set_yscale('log')\n", + "\n", + "ax1.grid()\n", + "ax2.grid()\n", + "\n", + "plt.tight_layout()\n", + "plt.savefig(\"HF_IA_vs_E_plot.png\", dpi=fig.dpi)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Search for superallowed decays in $\\epsilon/\\beta^{+}$ decay" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Superallowed 0+->0+ decays in EC/BP nuclides\n", + "\n", + "superallowed = []\n", + "for k,v, in ecbp_pairs.items():\n", + " decay = e.get_ecbp(edata, str(k[0]), int(k[3]), units='best')\n", + " for k_decay, v_decay in decay.items():\n", + " # Firm J and pi\n", + " if int(k_decay[9])==1 and int(k_decay[11])==1:\n", + " # Find parents with Jpi=0+\n", + " if int(2*k_decay[8])==0 and int(k_decay[10])==1:\n", + " # Find daughter with Jpi=0+\n", + " for daughter in v_decay:\n", + " # Search for ground states only\n", + " #if int(daughter[0])==0:\n", + " # Firm J and pi:\n", + " if int(daughter[2]==1) and int(daughter[4])==1 and int(daughter[6])==1:\n", + " # Find J=0 and pi=+:\n", + " try:\n", + " if int(2*daughter[3])==0 and int(daughter[5])==1 and daughter[17]<4.0:\n", + " print(k[0], k[1], k[2], k[3], k[4], k_decay[3], k_decay[4], \\\n", + " k_decay[5], daughter[0], daughter[1], daughter[3], \\\n", + " daughter[5], daughter[17], daughter[18], daughter[21], \\\n", + " daughter[22], daughter[23])\n", + " superallowed.append([k[0], k[1], k[2], k[3], k[4], k_decay[3], k_decay[4], \\\n", + " k_decay[5], daughter[0], daughter[1], daughter[3], \\\n", + " daughter[5], daughter[17], daughter[18], daughter[21], \\\n", + " daughter[22], daughter[23]])\n", + " except TypeError:\n", + " pass" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for l in superallowed: print(l)\n", + "print(len(superallowed))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sa = [i for i in superallowed if i[12] is not None]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(len(sa))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sa_array = np.array(sa)\n", + "A = sa_array[:,2].astype(int)\n", + "logft = sa_array[:,12].astype(float)\n", + "d_logft = sa_array[:,13].astype(float)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "fig, ax = plt.subplots(figsize=(9,6))\n", + "\n", + "ax.errorbar(A, logft, yerr=d_logft, color='k', fmt='o', capsize=5, label=r'Superallowed $0^{+} \\rightarrow 0^{+}$')\n", + "#ax.scatter(A, logft, color='k', label=r'Superallowed $0^{+} \\rightarrow 0^{+}$')\n", + "ax.legend(loc='best', fontsize=15)\n", + "ax.set_xlabel(r'$A$', size=20)\n", + "ax.set_ylabel(r'log$(ft)$', size=20)\n", + "\n", + "plt.grid()\n", + "plt.tight_layout()\n", + "plt.savefig(\"superallowed.png\", dpi=fig.dpi)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "x=e.get_ecbp(edata, \"S28\", 0, units=\"best\")\n", + "for k,v in x.items():\n", + " print(k,v)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Search for first-forbidden transitions in $\\beta^{-}$ decay" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Find 1st forbidden transitions based on a.m. selection rules and compare logft values\n", + "forbidden_first_list = []\n", + "forbidden_first=e.find_forbidden(edata,'1F',mode='ECBP')\n", + "for k,v in forbidden_first.items():\n", + " for vv in v:\n", + " if vv[7] is not None:\n", + " forbidden_first_list.append([k[11], vv[7], vv[8]])\n", + "\n", + "forbidden_firstu_list = []\n", + "forbidden_firstu=e.find_forbidden(edata,'1UF',mode='ECBP')\n", + "for k,v in forbidden_firstu.items():\n", + " for vv in v:\n", + " if vv[7] is not None:\n", + " forbidden_firstu_list.append([k[11], vv[7], vv[8]])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "first_array = np.array(forbidden_first_list)\n", + "f_A = first_array[:,0].astype(int)\n", + "f_logft = first_array[:,1].astype(float)\n", + "f_dlogft = first_array[:,2].astype(float)\n", + "\n", + "firstu_array = np.array(forbidden_firstu_list)\n", + "fu_A = firstu_array[:,0].astype(int)\n", + "fu_logft = firstu_array[:,1].astype(float)\n", + "fu_dlogft = firstu_array[:,2].astype(float)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "fig, ax = plt.subplots(figsize=(9,6))\n", + "\n", + "ax.errorbar(f_A, f_logft, yerr=f_dlogft, color='k', fmt='o', capsize=5, alpha=0.5, label=r'First forbidden')\n", + "ax.errorbar(fu_A, fu_logft, yerr=fu_dlogft, color='r', fmt='o', capsize=5, alpha=0.5, label=r'First forbidden unique')\n", + "ax.legend(loc='best', fontsize=15)\n", + "ax.set_xlabel(r'$A$', size=20)\n", + "ax.set_ylabel(r'log$(ft)$', size=20)\n", + "\n", + "ax.hlines(y=6, xmin=0, xmax=260, linewidth=2, color='b')\n", + "ax.hlines(y=9, xmin=0, xmax=260, linewidth=2, color='b')\n", + "\n", + "plt.grid()\n", + "plt.tight_layout()\n", + "plt.savefig(\"forbidden.png\", dpi=fig.dpi)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/pyproject.toml b/pyproject.toml index 6bb2acc..c41629a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,6 +41,7 @@ where = ["."] [tool.mypy] python_version = "3.10" strict = true +exclude = ["notebooks", "scripts/notebook_plots.py"] [tool.black] line-length = 88