beta_shape_pnpi/notebooks/exp_calibration.ipynb

Experimental Data Calibration Pipeline

Interactive workflow for processing raw CAMAC MCA binary spectra: read -> subtract -> calibrate -> livetime -> validate -> export

Prerequisites

• Raw binary spectrum files (A*.DAT format)
• Am-241 calibration file (CAM1.DAT)
• Pulser frequency and real time (from calibration run metadata)

In [ ]:

import matplotlib
# Use non-interactive backend only when no display is available (e.g. jupyter execute)
import matplotlib as mpl
try:
    import matplotlib.pyplot as plt
    # Check if we're in a headless environment
    import os
    if not os.environ.get("DISPLAY") and not os.environ.get("WAYLAND_DISPLAY"):
        mpl.use("Agg")
except ImportError:
    mpl.use("Agg")
    import matplotlib.pyplot as plt
import numpy as np

from exp_data.pipeline import process_spectrum, process_run_sequence
from exp_data.root_io import write_spectrum, read_spectrum
from exp_data.livetime import LivetimeDetermination

Configuration

In [ ]:

# Paths
RAW_FILE = "../data/raw/A1"
CALIBRATION_FILE = "../data/raw/CAM1.DAT"
RUN_DIR = "../data/raw/"

# Pulser parameters (from calibration run metadata)
PULSER_FREQUENCY_HZ = 100.0
REAL_TIME_SEC = 3600.0

Raw Spectrum Visualization

Visualize raw spectra before calibration to identify peaks and verify calibration source.

In [ ]:

# Visualize raw spectra before calibration
# This helps identify what calibration source is actually present

# Detect header size from file
def detect_header_size(filepath):
    """Detect MCA header size by checking for binary pattern."""
    with open(filepath, 'rb') as f:
        header = f.read(128)
    # Check if first 32 bytes look like a header (small ints)
    # and bytes 32-128 also look like header
    # Simple heuristic: if file is exactly 32 + 4096*4 = 16412, header=32
    # If file is 128 + 4096*4 = 16512, header=128
    import os
    size = os.path.getsize(filepath)
    expected_32 = 32 + 4096 * 4
    expected_128 = 128 + 4096 * 4
    if size == expected_32:
        return 32
    elif size == expected_128:
        return 128
    else:
        # Try to detect by checking if byte 32 looks like an int
        with open(filepath, 'rb') as f:
            f.seek(32)
            test = f.read(4)
            return 32  # default

# --- Raw run spectrum ---
header_size = detect_header_size(RAW_FILE)
with open(RAW_FILE, 'rb') as f:
    f.read(header_size)
    raw_counts = np.fromfile(f, dtype=np.int32, count=4096)

fig, axes = plt.subplots(2, 2, figsize=(14, 8))

ax = axes[0, 0]
ax.plot(raw_counts, 'b-', linewidth=0.5)
ax.set_title('Raw Run Spectrum (counts vs channel)')
ax.set_xlabel('Channel')
ax.set_ylabel('Counts')
ax.set_yscale('log')
ax.set_xlim(0, 600)

ax = axes[0, 1]
ax.plot(raw_counts, 'b-', linewidth=0.5)
ax.set_title('Raw Run Spectrum (full range)')
ax.set_xlabel('Channel')
ax.set_ylabel('Counts')
ax.set_yscale('log')

# --- Calibration spectrum ---
with open(CALIBRATION_FILE, 'rb') as f:
    f.read(128)  # CAM files have 128-byte header
    cal_counts = np.fromfile(f, dtype=np.int32, count=4096)

ax = axes[1, 0]
ax.plot(cal_counts, 'r-', linewidth=0.5)
ax.set_title('Calibration Spectrum (counts vs channel)')
ax.set_xlabel('Channel')
ax.set_ylabel('Counts')
ax.set_xlim(0, 600)
ax.set_yscale('log')

ax = axes[1, 1]
ax.plot(cal_counts, 'r-', linewidth=0.5)
ax.set_title('Calibration Spectrum (full range)')
ax.set_xlabel('Channel')
ax.set_ylabel('Counts')
ax.set_yscale('log')

plt.tight_layout()
plt.show()

# Print some stats
print(f'Raw run: {len(raw_counts)} channels, max={raw_counts.max()} at ch {np.argmax(raw_counts)}')
print(f'Calibration: {len(cal_counts)} channels, max={cal_counts.max()} at ch {np.argmax(cal_counts)}')

Process Single Spectrum

In [ ]:

result = process_spectrum(
    raw_path=RAW_FILE,
    calibration_path=CALIBRATION_FILE,
    pulser_frequency_hz=PULSER_FREQUENCY_HZ,
    real_time_sec=REAL_TIME_SEC,
)

print(f"Validation: {'PASSED' if result.validation_passed else 'FAILED'}")
print(f"Calibration chi2/dof: {result.calibration.chi2_per_dof:.4f}")
for msg in result.validation_messages:
    print(f"  {msg}")
if result.livetime:
    print(f"Live-time fraction: {result.livetime.live_time_fraction:.4f}")
    print(f"Pulser peak: {result.livetime.pulser_peak_center_keV:.1f} keV")

Plot Calibrated Spectrum

In [ ]:

plt.figure(figsize=(10, 6))
plt.plot(result.spectrum.energies, result.spectrum.counts, 'b-', label='Counts')
plt.xlabel('Energy (keV)')
plt.ylabel('Counts')
plt.yscale('log')
plt.title('Calibrated Beta Spectrum')
plt.legend()
plt.tight_layout()
plt.show()

Process Run Sequence

In [ ]:

results = process_run_sequence(
    run_dir=RUN_DIR,
    calibration_path=CALIBRATION_FILE,
    prefix="A",
    pulser_frequency_hz=PULSER_FREQUENCY_HZ,
    real_time_sec=REAL_TIME_SEC,
)

print(f"Processed {len(results)} runs")
for i, r in enumerate(results):
    print(f"  Run {i+1}: {'PASSED' if r.validation_passed else 'FAILED'} "
          f"(chi2/dof={r.calibration.chi2_per_dof:.4f})")

Export to ROOT

In [ ]:

write_spectrum("output.root", result.spectrum)
print("Exported to output.root")

Read Back from ROOT

In [ ]:

read_back = read_spectrum("output.root")
print(f"Read {len(read_back.energies)} channels")
print(f"Source: {read_back.source}")
print(f"Run ID: {read_back.run_id}")
print(f"Calibration: {read_back.metadata.get('calibration', {})}")

# Verify round-trip
assert np.allclose(read_back.energies, result.spectrum.energies)
assert np.allclose(read_back.counts, result.spectrum.counts)
print("Round-trip verified")

Livetime Determination

In [ ]:

# Direct livetime determination from calibration file
livetime_result = LivetimeDetermination.from_file(
    filepath=CALIBRATION_FILE,
    pulser_frequency_hz=PULSER_FREQUENCY_HZ,
    real_time_sec=REAL_TIME_SEC,
    source="pulser",
)

print(f"Live-time fraction: {livetime_result.live_time_fraction:.4f}")
print(f"Observed pulser rate: {livetime_result.observed_pulser_rate:.1f} Hz")
print(f"Pulser peak: {livetime_result.pulser_peak_center_keV:.1f} keV "
      f"(sigma={livetime_result.pulser_peak_sigma_keV:.2f} keV, "
      f"chi2/dof={livetime_result.pulser_peak_chi2:.2f})")