fix(scripts): list min/max calls, typed mains, drop unused header reads
This commit is contained in:
@@ -94,8 +94,6 @@
|
|||||||
" endpoint_MeV=0.294,\n",
|
" endpoint_MeV=0.294,\n",
|
||||||
" transition_type=\"F2\",\n",
|
" transition_type=\"F2\",\n",
|
||||||
" e_step_MeV=0.001,\n",
|
" e_step_MeV=0.001,\n",
|
||||||
" use_phase_space=True,\n",
|
|
||||||
" use_fermi=True,\n",
|
|
||||||
" use_screening=True,\n",
|
" use_screening=True,\n",
|
||||||
" use_finite_size=True,\n",
|
" use_finite_size=True,\n",
|
||||||
" use_charge_dist=True,\n",
|
" use_charge_dist=True,\n",
|
||||||
@@ -410,13 +408,16 @@
|
|||||||
}
|
}
|
||||||
],
|
],
|
||||||
"source": [
|
"source": [
|
||||||
|
"multi_analyzer = BetaSpectrumAnalyzer(multi_spectrum, multi_config)\n",
|
||||||
|
"\n",
|
||||||
"fig, ax = plt.subplots(figsize=(12, 6))\n",
|
"fig, ax = plt.subplots(figsize=(12, 6))\n",
|
||||||
"ax.plot(T_keV, total, \"k-\", linewidth=2.5, label=\"Total spectrum\")\n",
|
"ax.plot(T_keV, total, \"k-\", linewidth=2.5, label=\"Total spectrum\")\n",
|
||||||
"\n",
|
"\n",
|
||||||
"branch_spectra = multi_spectrum.get_branch_normalized_spectra(W)\n",
|
"# Physical branch spectra: core (phase_space x fermi) x all enabled corrections\n",
|
||||||
"for i, bs in enumerate(branch_spectra):\n",
|
"for i, b in enumerate(branches):\n",
|
||||||
" ax.plot(T_keV, bs, alpha=0.6, linewidth=1.5,\n",
|
" physical = multi_analyzer.physical_branch_spectrum(i)\n",
|
||||||
" label=\"Branch {} (E0={:.0f} keV)\".format(i+1, branches[i].endpoint_MeV*1000))\n",
|
" ax.plot(T_keV, physical, alpha=0.6, linewidth=1.5,\n",
|
||||||
|
" label=\"Branch {} (E0={:.0f} keV)\".format(i+1, b.endpoint_MeV*1000))\n",
|
||||||
"\n",
|
"\n",
|
||||||
"for i, b in enumerate(branches):\n",
|
"for i, b in enumerate(branches):\n",
|
||||||
" br_end_keV = W_to_T(T_to_W(b.endpoint_MeV))\n",
|
" br_end_keV = W_to_T(T_to_W(b.endpoint_MeV))\n",
|
||||||
@@ -430,7 +431,7 @@
|
|||||||
"ax.set_xlim(0, W0_max_keV)\n",
|
"ax.set_xlim(0, W0_max_keV)\n",
|
||||||
"ax.grid(True, alpha=0.3)\n",
|
"ax.grid(True, alpha=0.3)\n",
|
||||||
"plt.tight_layout()\n",
|
"plt.tight_layout()\n",
|
||||||
"plt.show()"
|
"plt.show()\n"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
@@ -541,8 +542,8 @@
|
|||||||
}
|
}
|
||||||
],
|
],
|
||||||
"source": [
|
"source": [
|
||||||
"analyzer = BetaSpectrumAnalyzer(multi_spectrum, multi_config)\n",
|
"analyzer = multi_analyzer\n",
|
||||||
"analyzer.plot_analysis(show_components=True)"
|
"analyzer.plot_analysis(show_components=True)\n"
|
||||||
]
|
]
|
||||||
},
|
},
|
||||||
{
|
{
|
||||||
|
|||||||
@@ -81,11 +81,7 @@
|
|||||||
"\n",
|
"\n",
|
||||||
"# Detect header size from file\n",
|
"# Detect header size from file\n",
|
||||||
"def detect_header_size(filepath):\n",
|
"def detect_header_size(filepath):\n",
|
||||||
" \"\"\"Detect MCA header size by checking for binary pattern.\"\"\"\n",
|
" \"\"\"Detect MCA header size from the file size.\"\"\"\n",
|
||||||
" with open(filepath, 'rb') as f:\n",
|
|
||||||
" header = f.read(128)\n",
|
|
||||||
" # Check if first 32 bytes look like a header (small ints)\n",
|
|
||||||
" # and bytes 32-128 also look like header\n",
|
|
||||||
" # Simple heuristic: if file is exactly 32 + 4096*4 = 16412, header=32\n",
|
" # Simple heuristic: if file is exactly 32 + 4096*4 = 16412, header=32\n",
|
||||||
" # If file is 128 + 4096*4 = 16512, header=128\n",
|
" # If file is 128 + 4096*4 = 16512, header=128\n",
|
||||||
" import os\n",
|
" import os\n",
|
||||||
@@ -97,10 +93,6 @@
|
|||||||
" elif size == expected_128:\n",
|
" elif size == expected_128:\n",
|
||||||
" return 128\n",
|
" return 128\n",
|
||||||
" else:\n",
|
" else:\n",
|
||||||
" # Try to detect by checking if byte 32 looks like an int\n",
|
|
||||||
" with open(filepath, 'rb') as f:\n",
|
|
||||||
" f.seek(32)\n",
|
|
||||||
" test = f.read(4)\n",
|
|
||||||
" return 32 # default\n",
|
" return 32 # default\n",
|
||||||
"\n",
|
"\n",
|
||||||
"# --- Raw run spectrum ---\n",
|
"# --- Raw run spectrum ---\n",
|
||||||
|
|||||||
+36
-15
@@ -17,6 +17,7 @@ from __future__ import annotations
|
|||||||
|
|
||||||
import sys
|
import sys
|
||||||
from pathlib import Path
|
from pathlib import Path
|
||||||
|
from typing import cast
|
||||||
|
|
||||||
import numpy as np
|
import numpy as np
|
||||||
from scipy.stats import linregress
|
from scipy.stats import linregress
|
||||||
@@ -101,7 +102,10 @@ def estimate_endpoint_from_counts(
|
|||||||
|
|
||||||
# Linear fit: K = a + b * E, endpoint = -a/b
|
# Linear fit: K = a + b * E, endpoint = -a/b
|
||||||
try:
|
try:
|
||||||
slope, intercept, r_val, _, _ = linregress(E_sel, K)[:5]
|
lin = cast("tuple[float, float, float, float, float]", linregress(E_sel, K))
|
||||||
|
slope = lin[0]
|
||||||
|
intercept = lin[1]
|
||||||
|
r_val = lin[2]
|
||||||
if slope < 0 and abs(r_val) > 0.3:
|
if slope < 0 and abs(r_val) > 0.3:
|
||||||
endpoint = -intercept / slope
|
endpoint = -intercept / slope
|
||||||
# Clamp to reasonable range (Q ± 100 keV)
|
# Clamp to reasonable range (Q ± 100 keV)
|
||||||
@@ -129,7 +133,7 @@ def compute_regional_rates(counts: np.ndarray, energies_keV: np.ndarray) -> dict
|
|||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
def main():
|
def main() -> None:
|
||||||
"""Main analysis pipeline."""
|
"""Main analysis pipeline."""
|
||||||
data_dir = Path(__file__).parent.parent / "data" / "raw"
|
data_dir = Path(__file__).parent.parent / "data" / "raw"
|
||||||
output_dir = Path(__file__).parent.parent / "output"
|
output_dir = Path(__file__).parent.parent / "output"
|
||||||
@@ -168,8 +172,6 @@ def main():
|
|||||||
print("\nSubtracting consecutive runs...")
|
print("\nSubtracting consecutive runs...")
|
||||||
individual_runs = []
|
individual_runs = []
|
||||||
for i in range(1, len(cumulative_specs)):
|
for i in range(1, len(cumulative_specs)):
|
||||||
parent_a = cumulative_specs[i].run_id
|
|
||||||
parent_b = cumulative_specs[i - 1].run_id
|
|
||||||
diff = subtract_spectra(
|
diff = subtract_spectra(
|
||||||
cumulative_a=cumulative_specs[i],
|
cumulative_a=cumulative_specs[i],
|
||||||
cumulative_b=cumulative_specs[i - 1],
|
cumulative_b=cumulative_specs[i - 1],
|
||||||
@@ -215,33 +217,52 @@ def main():
|
|||||||
print(f"\n{'=' * 60}")
|
print(f"\n{'=' * 60}")
|
||||||
print(f"Analysis of {len(run_ids)} individual Tc-99 runs")
|
print(f"Analysis of {len(run_ids)} individual Tc-99 runs")
|
||||||
print(f"{'=' * 60}")
|
print(f"{'=' * 60}")
|
||||||
print(f"Run range: {run_ids.min()} - {run_ids.max()}")
|
print(f"Run range: {min(run_ids)} - {max(run_ids)}")
|
||||||
print(f"Total counts: {total_rates.min():.0f} - {total_rates.max():.0f}")
|
print(f"Total counts: {min(total_rates):.0f} - {max(total_rates):.0f}")
|
||||||
print(f"Low-energy fraction: {low_fractions.min():.3f} - {low_fractions.max():.3f}")
|
print(
|
||||||
print(f"Endpoint estimate: {endpoints.min():.1f} - {endpoints.max():.1f} keV")
|
f"Low-energy fraction: {min(low_fractions):.3f} " f"- {max(low_fractions):.3f}"
|
||||||
|
)
|
||||||
|
print(f"Endpoint estimate: {min(endpoints):.1f} - {max(endpoints):.1f} keV")
|
||||||
|
|
||||||
# Endpoint drift
|
# Endpoint drift
|
||||||
if len(run_ids) > 10:
|
if len(run_ids) > 10:
|
||||||
slope, intercept, r_val, _, _ = linregress(run_ids, endpoints)
|
lin_e = cast(
|
||||||
print(f"\nEndpoint drift: {slope:.3f} keV/run (R²={r_val**2:.4f})")
|
"tuple[float, float, float, float, float]",
|
||||||
print(f" Start (run 1): {intercept + slope * run_ids[0]:.1f} keV")
|
linregress(run_ids, endpoints),
|
||||||
print(f" End (run {run_ids[-1]}): {intercept + slope * run_ids[-1]:.1f} keV")
|
)
|
||||||
|
slope_e = lin_e[0]
|
||||||
|
intercept_e = lin_e[1]
|
||||||
|
r_e = lin_e[2]
|
||||||
|
print(f"\nEndpoint drift: {slope_e:.3f} keV/run (R²={r_e ** 2:.4f})")
|
||||||
|
print(f" Start (run 1): {intercept_e + slope_e * run_ids[0]:.1f} keV")
|
||||||
|
print(
|
||||||
|
f" End (run {run_ids[-1]}): {intercept_e + slope_e * run_ids[-1]:.1f} keV"
|
||||||
|
)
|
||||||
|
|
||||||
# Rate drift
|
# Rate drift
|
||||||
slope, intercept, r_val, _, _ = linregress(run_ids, total_rates)
|
lin_r = cast(
|
||||||
print(f"\nTotal count rate drift: {slope:.0f} counts/run (R²={r_val**2:.4f})")
|
"tuple[float, float, float, float, float]", linregress(run_ids, total_rates)
|
||||||
|
)
|
||||||
|
slope_r = lin_r[0]
|
||||||
|
r_r = lin_r[2]
|
||||||
|
print(f"\nTotal count rate drift: {slope_r:.0f} counts/run (R²={r_r ** 2:.4f})")
|
||||||
|
|
||||||
# Save results
|
# Save results
|
||||||
results_file = output_dir / "rate_analysis.csv"
|
results_file = output_dir / "rate_analysis.csv"
|
||||||
|
try:
|
||||||
with open(results_file, "w") as f:
|
with open(results_file, "w") as f:
|
||||||
f.write(
|
f.write(
|
||||||
"run,total_counts,low_energy_counts,high_energy_counts,low_fraction,endpoint_keV\n"
|
"run,total_counts,low_energy_counts,high_energy_counts,"
|
||||||
|
"low_fraction,endpoint_keV\n"
|
||||||
)
|
)
|
||||||
for i in range(len(run_ids)):
|
for i in range(len(run_ids)):
|
||||||
f.write(
|
f.write(
|
||||||
f"{run_ids[i]},{total_rates[i]:.0f},{low_rates[i]:.0f},"
|
f"{run_ids[i]},{total_rates[i]:.0f},{low_rates[i]:.0f},"
|
||||||
f"{high_rates[i]:.0f},{low_fractions[i]:.4f},{endpoints[i]:.1f}\n"
|
f"{high_rates[i]:.0f},{low_fractions[i]:.4f},{endpoints[i]:.1f}\n"
|
||||||
)
|
)
|
||||||
|
except OSError as exc:
|
||||||
|
print(f"Cannot write {results_file}: {exc}")
|
||||||
|
raise SystemExit(1) from exc
|
||||||
print(f"\nResults saved to {results_file}")
|
print(f"\nResults saved to {results_file}")
|
||||||
|
|
||||||
# Plot
|
# Plot
|
||||||
|
|||||||
@@ -471,7 +471,7 @@ CHECKS = [
|
|||||||
]
|
]
|
||||||
|
|
||||||
|
|
||||||
def main():
|
def main() -> int:
|
||||||
parser = argparse.ArgumentParser(description="Check documentation consistency")
|
parser = argparse.ArgumentParser(description="Check documentation consistency")
|
||||||
parser.add_argument(
|
parser.add_argument(
|
||||||
"--strict", action="store_true", help="Treat warnings as errors (exit 1)"
|
"--strict", action="store_true", help="Treat warnings as errors (exit 1)"
|
||||||
|
|||||||
Reference in New Issue
Block a user