Chapter 3
Alpha/Beta Spectra and Energy Quantisation

Flow Algebra

flow "AlphaBeta_Spectra_EnergyQuantisation"

# ---------- STATES ----------
state:
  # nuclide / energetics
  nucleus.Z: int
  nucleus.N: int
  Q.value: float                     # decay Q for active branch (MeV)
  level.set: set[level]              # {E_ex, Jπ, branch}
  branch.table: map[channel, float]  # branching ratios (sum ≤ 1)
  channel.active: str?               # "alpha" | "beta_minus" | "beta_plus" | "EC" | "gamma"

  # sampled event kinematics (per tick if simulating)
  Ee: float                          # electron/positron kinetic energy
  Enu: float                         # neutrino energy (β channels)
  Eα: float                          # alpha kinetic energy
  Eγ: float                          # gamma energy
  Erecoil: float                     # nuclear recoil
  Elevel: float                      # selected excited level energy

  # spectral accumulators (ideal pre-detector)
  spec.alpha_lines: list[(E:float, weight:float)]
  spec.beta_e_hist: histogram        # electron/positron energy spectrum
  spec.gamma_lines: list[(E:float, weight:float)]
  spec.XR_Auger_lines: list[(E:float, weight:float)]  # from EC

  # detector model
  det.model: "gaussian" | "voigt"
  det.FWHM_keV: float                # resolution at reference energy
  det.ref_E_keV: float               # energy where FWHM_keV applies
  det.resolution_law: "sqrt" | "linear" | "constant"
  det.eta_tail: float                # low-energy tail fraction (0..1), optional

  # quantisation / calibration
  bin.Emin_keV: float
  bin.Emax_keV: float
  bin.N: int
  bin.edges: array[float]
  bin.centres: array[float]
  calib.gain: float                  # keV/channel
  calib.offset: float                # keV
  calib.nonlinearity: float          # small quadratic term (keV^-1), optional

  # physics modifiers
  Fermi.Zeff: float                  # effective Z for F(Z,E)
  shape_factor.kind: "allowed" | "first_forbidden" | "custom"
  shape_factor.C: fn(E)->float       # multiplicative spectral shape
  screening.S: fn(E)->float          # atomic screening correction
  radiative.R: fn(E)->float          # radiative (inner brems) correction
  endpoint.correction: fn()->float   # tiny shifts (recoil, screening) at endpoint

  # normalisation / runtime
  activity.Bq: float                 # events per second (for cps)
  live_time.s: float
  dead_time.frac: float              # 0..1
  pileup.rate: float                 # events/s (optional)

  # outputs (post-detector, quantised)
  spectrum.alpha_counts: array[float]  # binned counts
  spectrum.beta_counts: array[float]
  spectrum.gamma_counts: array[float]
  spectrum.total_counts: array[float]

# ---------- PARAMETERS ----------
params:
  me_keV: 510.99895                  # electron rest mass
  alpha_mass_keV: 3727.379           # alpha rest mass
  c2: 1.0                            # MeV units (natural)
  N_MC: int                          # Monte Carlo samples for spectra build
  voigt.eta: float                   # Lorentzian fraction if det.model="voigt"

# ---------- INVARIANTS ----------
invariant:
  - sum(branch.table.values) ≤ 1 + 1e-9
  - for α:  Q.value = Eα + Erecoil + Elevel    # (per event)
  - for β−: Q.value = Ee + Enu + Erecoil + Elevel
  - for β+: Q.value = Ee + Enu + 2*params.me_keV/1000 + Erecoil + Elevel
  - Ee ≥ 0, Eα ≥ 0, Eγ ≥ 0, Erecoil ≥ 0, Elevel ≥ 0
  - β endpoints: E0 = Q.value - Elevel - endpoint.correction() ≥ 0
  - histogram area ∝ expected decays after live time and dead-time losses

# ---------- TRANSFORMS ----------
transform:

  initialise_bins:
    do:
      bin.edges   = linspace(bin.Emin_keV, bin.Emax_keV, bin.N+1)
      bin.centres = 0.5*(bin.edges[:-1] + bin.edges[1:])

  choose_channel:
    do:
      channel.active = random_weighted_pick(branch.table)
      if channel.active == null: channel.active = "gamma"  # if pure isomeric γ

  pick_level:
    do:
      # sample nuclear level from level.set, weighted by branch × level.branch
      L = random_weighted_pick({lvl: lvl.branch for lvl in level.set})
      Elevel = L.E_ex

  alpha_event:
    pre:
      require channel.active == "alpha"
    do:
      # two-body kinematics → sharp line, split with recoil
      Eα = Q.value - Elevel - Erecoil_alpha(Q.value, nucleus.Z, nucleus.N)
      append(spec.alpha_lines, (E=to_keV(Eα), weight=1.0))

  beta_event:
    pre:
      require channel.active in ["beta_minus","beta_plus"]
    do:
      E0 = Q.value - Elevel - endpoint.correction()
      require E0 > 0
      # sample Ee from allowed/forbidden spectrum
      # S(E) ∝ F(Z,E) * p * (E+me)^k * (E0 - E)^2 * shape_factor.C(E) * screening.S(E) * radiative.R(E)
      Ee = sample_beta_energy(E0, Fermi.Zeff, shape_factor, screening, radiative, kind=channel.active)
      add_to_hist(spec.beta_e_hist, to_keV(Ee), weight=1.0)

  EC_event:   # electron capture → characteristic X-ray/Auger lines
    pre:
      require channel.active == "EC"
    do:
      lines = atomic_relaxation_lines(nucleus.Z)  # Kα, Kβ, L..., Auger
      for (E, w) in lines:
        append(spec.XR_Auger_lines, (E=E, weight=w))

  gamma_event:
    pre:
      require channel.active == "gamma"
    do:
      Eγ = Elevel   # de-excitation
      append(spec.gamma_lines, (E=to_keV(Eγ), weight=1.0))

  build_spectra_MC:
    do:
      for i in range(params.N_MC):
        choose_channel()
        pick_level()
        match channel.active:
          - "alpha":        alpha_event()
          - "beta_minus":   beta_event()
          - "beta_plus":    beta_event()
          - "EC":           EC_event()
          - "gamma":        gamma_event()

  apply_detector_response:
    do:
      # 1) start with zeros
      spectrum.alpha_counts = zeros(bin.N)
      spectrum.beta_counts  = zeros(bin.N)
      spectrum.gamma_counts = zeros(bin.N)

      # 2) convolve ideal lines/spectra with detector kernel (energy-dependent σ)
      for (E, w) in spec.alpha_lines:
        K = detector_kernel(E, bin.centres, det.FWHM_keV, det.ref_E_keV, det.resolution_law, det.model, params.voigt.eta, det.eta_tail)
        spectrum.alpha_counts += w * K

      # β: take histogram then smooth by kernel per bin centre
      spectrum.beta_counts += convolve_hist(spec.beta_e_hist, bin.centres, det.FWHM_keV, det.ref_E_keV, det.resolution_law, det.model, params.voigt.eta, det.eta_tail)

      # γ & EC lines
      for (E, w) in spec.gamma_lines:
        K = detector_kernel(E, bin.centres, det.FWHM_keV, det.ref_E_keV, det.resolution_law, det.model, params.voigt.eta, det.eta_tail)
        spectrum.gamma_counts += w * K
      for (E, w) in spec.XR_Auger_lines:
        K = detector_kernel(E, bin.centres, det.FWHM_keV, det.ref_E_keV, det.resolution_law, det.model, params.voigt.eta, det.eta_tail)
        spectrum.gamma_counts += w * K   # included in same panel by default

  calibrate_and_quantise:
    do:
      # Map true keV → channel space, then back via gain/offset/nonlinearity to mimic ADC quantisation
      # Here we assume bins are already in keV; nonlinearity can warp counts slightly across bins if desired.
      if calib.nonlinearity != 0:
        warp_bins(bin.centres, calib.nonlinearity)

  scale_counts_to_cps:
    do:
      scale = activity.Bq * live_time.s * (1 - dead_time.frac) / params.N_MC
      spectrum.alpha_counts *= scale
      spectrum.beta_counts  *= scale
      spectrum.gamma_counts *= scale
      spectrum.total_counts  = spectrum.alpha_counts + spectrum.beta_counts + spectrum.gamma_counts

# ---------- MONITORS ----------
monitor:
  alpha_line_index:
    observe: [spec.alpha_lines]
    report: "Expected α line positions (keV) and relative weights before detector broadening."

  beta_endpoint_check:
    observe: [Q.value, Elevel]
    report: "Computed β endpoint E0 per level; ensure histograms cut off at E0 (post-kernel may bleed slightly)."

  resolution_health:
    observe: [det.FWHM_keV, det.resolution_law]
    report: "Check σ(E) scaling; too broad smears discrete lines, too narrow hides realism."

  normalisation:
    observe: [activity.Bq, live_time.s, dead_time.frac]
    report: "Counts should scale linearly with activity×live_time×(1-dead_time)."

# ---------- NOTES ----------
note:
  - α spectra are **discrete**: two-body kinematics → sharp lines; detector turns them into Gaussians/Voigts.
  - β spectra are **continuous** up to endpoint E0; shape ∝ Fermi function × phase space × (forbiddenness) × small corrections.
  - EC yields **characteristic X-ray/Auger** lines; include in γ panel or separate.
  - Recoil and excited-level energy subtract from particle kinetic energy, shifting α lines and β endpoints.
  - Detector resolution often follows σ(E) ∝ √E; choose det.resolution_law="sqrt" for scintillators/semiconductors.
  - Apply small endpoint corrections (screening, radiative, recoil) if comparing to high-precision data.
  - For β+ include the 2me term in Q accounting; if modelling annihilation, add 511 keV γ lines (2× or 3× topology).