CMA Spectral Unfolder

Multi-species spectral unfolding from scintillator stack detector data using CMA-ES (Covariance Matrix Adaptation Evolution Strategy), with Huber-robust fidelity, Tikhonov / logistic smoothness regularisation, and optional per-detector calibration tuning.

This code reconstructs particle energy spectra from a stack of scintillator detectors by inverting a pre-computed FLUKA Monte Carlo response matrix (RM). The inverse problem is ill-posed: the code solves it by minimising a regularised objective function using the derivative-free CMA-ES optimizer, supporting up to three particle species simultaneously.

Physics Background

The Inverse Problem

A stack of N scintillator slabs records the energy deposited by a broad-spectrum particle beam. The measured signal in detector channel \(d_j\) (\(j = 1 \ldots N_\text{FLUKA}\)) is a linear superposition of mono-energetic responses:

where \(R_j(E_i)\) is the response matrix element (energy deposited in detector \(j\) by a particle of energy \(E_i\)), and \(S_i = (dN/dE)_i\) is the unknown spectrum evaluated at the \(i\)-th energy bin. In matrix form: \(\mathbf{d} = \mathbf{R}\,\mathbf{s}\).

Why CMA-ES?

The response matrix is typically rank-deficient or severely ill-conditioned (condition numbers of \(10^6\)–\(10^{12}\) are common for thick-stack calorimeters). Classical linear inversions (SVD truncation, Tikhonov) struggle with positivity constraints and multi-species superposition. CMA-ES is a gradient-free evolutionary strategy that handles:

• Non-convex, multi-modal landscapes
• Box constraints (positivity via log-space parameterisation)
• Heterogeneous variable types (spectra + calibration factors)
• Arbitrary regularisation terms without differentiability requirements

Multi-Species Unfolding

When the detector is sensitive to multiple particle types (e.g. photons + electrons), the measured signal is:

where \(k\) indexes the species, each with its own RM and energy grid. The optimiser solves for all species simultaneously, which breaks the degeneracy that arises from fitting a single RM to a mixed signal.

Architecture

The code is organised into six classes and a set of free functions:

Getting Started

# Install dependencies
pip install numpy pandas matplotlib pillow scipy cma

# Clone and enter the project
git clone https://github.com/your-org/cma-unfolder.git
cd cma-unfolder

cma-unfolder/
├── CMA_optimizer.py          # Main script
├── images/
│   └── raw_data.tiff         # Scintillator image
├── RM/
│   └── Response_matrix_double_population/
│       ├── Response_matrix_p/ # Photon RM files
│       │   ├── 0.05_Energy.txt
│       │   ├── 0.05_Spectrum.txt
│       │   ├── 0.05_FLUKA.txt
│       │   ├── 0.1_Energy.txt
│       │   └── ...
│       └── Response_matrix_e/ # Electron RM files
│           └── ...
├── Error/
│   ├── mean1_VAC.txt          # Error level 1 (signal > 10%)
│   ├── mean2_VAC.txt          # Error level 2 (1% < signal < 10%)
│   ├── mean3_VAC.txt          # Error level 3 (0.1% < signal < 1%)
│   └── mean4_VAC.txt          # Error level 4 (signal < 0.1%)
└── README.md

Configuration — Config

All user-tunable parameters live in the Config class as class-level attributes. Modify this class to adapt the code to a new experimental setup. No constructor is needed — all values are set at import time.

Forward Model

For a single species with spectrum \(\{S_i\}_{i=1}^{N_\text{guess}}\) in log₁₀ space, the predicted detector signal is:

where:

• \(R_{ij}\) = FLUKA_tot[nearest[i], j] — normalised RM element
• \(f_i\) = FLUKA_fact[nearest[i]] — peak normalization factor (restores absolute scale)
• \(c_j\) = facts[j] — optional per-detector calibration correction (\(\approx 1\))

For \(K\) species the total signal is the sum: \(\hat{d}_j = \sum_{k=1}^{K} \hat{d}_j^{(k)}\). This is implemented in _calc_FLUKA_species() and is fully vectorised over energy bins.

Loss Function

The total objective minimised by CMA-ES is:

Fidelity Term

Relative squared residuals with Huber robustification:

where \(H_\delta(r) = \frac{1}{2}\min(|r|,\delta)^2 + \delta\,(|r| - \min(|r|,\delta))\) transitions from quadratic to linear at threshold \(\delta\). This suppresses the influence of outlier detector channels (e.g. damaged scintillators or channels with poor RM coverage).

Smoothness Penalties

Simple (Tikhonov-like)

When adaptive_smoothing = False:

This is the squared \(\ell_2\) norm of the second finite difference — a discrete approximation to the curvature penalty \(\int |S''(E)|^2 dE\).

Logistic-Weighted

When adaptive_smoothing = True, the penalty uses first differences with bin-dependent weights:

where \(x = i/(N-1)\) is the normalised bin index and \(\sigma\) is the logistic function. This reduces regularisation strength at the high-energy tail where the RM sensitivity drops and the spectrum is expected to fall sharply, avoiding the artificial flattening that uniform smoothing produces near the cutoff.

Calibration Factor Tuning

When tune_calibration = True, \(N_\text{FLUKA}\) additional variables are appended to the CMA vector, representing per-detector multiplicative corrections \(c_j\). They are initialised near 1.0 and constrained by a smooth quartic barrier:

where \(\text{sp}(x) = \frac{1}{k}\log(1 + e^{kx})\) is the softplus function with sharpness \(k = 300\). This allows small deviations from unity while strongly penalising large excursions.

CMA-ES Internals

The optimizer uses pycma with the following default settings:

The objective function is dispatched once before the loop via _build_objective_fn(), which selects among four variants (base / smooth / calib / smooth+calib) based on the active flags. This eliminates per-evaluation branching overhead.

Optimisation Vector Layout

x = [ S_photon (N_guess_p) | S_electron (N_guess_e) | facts (N_FLUKA) ]
      ←── spectral vars ──────────────────────────→  ←── optional ──→

API Reference

SpeciesConfig

Dataclass. Container for one particle species' RM and optimiser settings.

Constructor Parameters

Methods

resolve()

Must be called once before optimisation. Maps E_guess onto the RM energy axis via nearest-neighbour lookup (injective), sets nearest_indices and N_guess.

MultiSpeciesOptimizer

Constructor Parameters

Methods

run_CMA() → OptimizationResult

Execute the CMA-ES evolution loop. Returns a structured OptimizationResult.

calc_FLUKA(result) → ndarray

Reconstruct the total detector signal from an OptimizationResult. Useful for computing post-hoc residuals.

mmin(x) → float

Internal objective function called by CMA-ES. Unpacks x and delegates to the pre-selected variant.

OptimizationResult

Dataclass returned by run_CMA().

get(label) → (E_guess, spectrum_linear)

Convenience accessor. Raises KeyError if the label is not found.

E, spec = result.get('photon')
plt.loglog(E, spec)

DataProcessor

import_RM() → (Spectrum_tot, FLUKA_tot, FLUKA_fact, E)

Reads all *_FLUKA.txt / *_Energy.txt file triplets from Config.folder_path. Applies Config.factor, peak-normalises each row, and returns arrays sorted by ascending energy.

read_image(norm_flag=True, plot_flag=False) → ndarray

Opens Config.image_path, extracts integrated pixel intensity per ROI. If norm_flag is True, each value is divided by the ROI pixel count.

find_nearest(array, values) → ndarray

Injective nearest-neighbour mapping: returns one unique index per value. Raises ValueError if any value falls outside the array range.

ErrorAnalysis

Constructor

ErrorAnalysis(E_error, error_files, ddv, window_size=200)

Loads four error-level files and builds a 2-D error matrix over (normalised signal, energy), smoothed with a moving-average kernel of width window_size.

get_errors(x_vals, data_vals) → ndarray

Bilinear interpolation on the error matrix. data_vals is peak-normalised internally. Returns NaN where interpolation fails (out-of-bounds).

plot_error_results(E_guess, signal, y_errors, std)

Produces a publication-ready log–log plot with error bars combining unfolding uncertainty, systematic error (std), and Poisson counting noise.

Plotter

plot_results(Exp_FLUKA, FLUKA_sim)

Overlays experimental and reconstructed detector signals (normalised).

plot_spectrum(E_guess, sim)

Log–log plot of the unfolded spectrum.

Helper Functions

Tutorial: Single-Species Unfolding

The simplest use case: unfold a photon spectrum from 17 scintillator channels.

from CMA_optimizer import (
    Config, DataProcessor, SpeciesConfig,
    MultiSpeciesOptimizer, Plotter
)

# 1. Load RM
_, FLUKA_tot, FLUKA_fact, E = DataProcessor.import_RM()

# 2. Read experimental data
Exp_FLUKA = DataProcessor.read_image(norm_flag=True)
Exp_max = Exp_FLUKA.max()
Exp_FLUKA /= Exp_max

# 3. Configure species
species = SpeciesConfig(
    label='photon',
    E=E,
    E_guess=Config.E_guess.copy(),
    FLUKA_tot=FLUKA_tot,
    FLUKA_fact=FLUKA_fact,
    smooth_factor=5e-6,
)

# 4. Run optimizer
opt = MultiSpeciesOptimizer(
    species_list=[species],
    Exp_FLUKA=Exp_FLUKA,
    smoothing=True,
)
result = opt.run_CMA()

# 5. Plot
E_out, spec = result.get('photon')
Plotter.plot_spectrum(E_out, spec * Exp_max)

Tutorial: Multi-Species (Photon + Electron)

# After loading photon RM as above...

# Load electron RM by temporarily swapping folder_path
_backup = Config.folder_path
Config.folder_path = Config.folder_path_e
_, FLUKA_tot_e, FLUKA_fact_e, E_e = DataProcessor.import_RM()
Config.folder_path = _backup

species_e = SpeciesConfig(
    label='electron',
    E=E_e,
    E_guess=Config.E_guess_e,
    FLUKA_tot=FLUKA_tot_e,
    FLUKA_fact=FLUKA_fact_e,
    smooth_factor=1.3e-9,
    smooth_transition_start=0.30,
    smooth_min_weight=0.10,
    smooth_steepness=10,
)

opt = MultiSpeciesOptimizer(
    species_list=[species_p, species_e],
    Exp_FLUKA=Exp_FLUKA,
    smoothing=True,
)
result = opt.run_CMA()

# Access each species independently
E_p, spec_p = result.get('photon')
E_e, spec_e = result.get('electron')

Tutorial: Synthetic Closure Test

The main() function demonstrates a closure test: a known analytical spectrum is forward-modelled through the RM to produce synthetic data, then recovered by the optimizer. This validates the pipeline end-to-end.

from CMA_optimizer import make_synthetic_exp

# Define a known exponential photon spectrum
E_temp = 10   # MeV temperature
known = (1 / Config.E_guess) * np.exp(-Config.E_guess / E_temp)

# Generate synthetic detector signal
Exp_FLUKA = make_synthetic_exp(FLUKA_tot, FLUKA_fact, E, Config.E_guess, known)
Exp_max = Exp_FLUKA.max()
Exp_FLUKA /= Exp_max

# Run optimizer and compare result.spectra_linear[0] against `known`

Tutorial: Calibration Factor Tuning

opt = MultiSpeciesOptimizer(
    species_list=[species_p],
    Exp_FLUKA=Exp_FLUKA,
    tune_calibration=True,
    facts_penalty_weight=1e-2,
    smoothing=True,
)
result = opt.run_CMA()

# Inspect optimised calibration factors
print(result.facts)  # should be close to 1.0

Hyperparameter Guide

Diagnostics & Debugging

Convergence Checks

• Inspect result.es_result.stop() to see which CMA termination criterion fired.
• Use cma.plot() on the logger for convergence curves (objective value, σ, axis ratios).
• If maxiter is reached without convergence, increase it or widen sigma0.

Residual Analysis

FLUKA_sim = opt.calc_FLUKA(result)
residual = (Exp_FLUKA.flatten() - FLUKA_sim / FLUKA_sim.max())
plt.bar(range(Config.N_FLUKA), residual)
plt.ylabel('Relative residual')
plt.xlabel('Detector channel')
plt.show()

Common Failure Modes

Performance & Scaling

Runtimes are approximate for a single-core desktop (3 GHz). CMA-ES is inherently sequential in its covariance update, but the objective function evaluations within each generation are independent and could be parallelised via cma.CMAEvolutionStrategy's eval_parallel callback. The forward model itself (_calc_FLUKA_species) is fully vectorised with NumPy.

Extending the Code

Adding a Third Species

Create a third SpeciesConfig with its own RM and energy grid, then pass all three to MultiSpeciesOptimizer. The code supports up to 3 species natively.

Custom Regularisation

To add a new regularisation term (e.g. entropy, L1 sparsity), define a function with signature f(spectrum: ndarray) → float and add it to the appropriate _objective_* method. Create a new objective variant if needed and register it in _build_objective_fn().

Alternative Fidelity Metrics

The Huber loss can be replaced by any robust estimator (Cauchy, Tukey biweight) by modifying _huber_loss(). The relative-residual normalisation \((d - \hat{d})^2 / d^2\) can also be changed to absolute residuals or Poisson log-likelihood depending on the noise model.

Known Issues

• Config as mutable class attributes: swapping Config.folder_path to load a second RM is a side-effect pattern. A future refactor should pass the folder path as an argument to import_RM().
• ErrorAnalysis double-smoothing: the error matrix is smoothed once inside _create_error_matrix() and again in __init__. This is inherited from the original code and produces slightly more diffuse error estimates than a single pass.
• Injective nearest-neighbour in find_nearest: if two E_guess values map to the same RM index, the second is bumped to the next available index. This can shift bins at the edges of the energy grid. Prefer well-separated E_guess values.
• No automatic normalisation consistency check between Exp_FLUKA and the RM scale. The user must ensure consistent units.

License

This software is released under the PolyForm Noncommercial License 1.0.0. You may use, modify, and distribute the code for any non-commercial purpose. Commercial use requires a separate license from the author.

© 2025 G. Fauvel — CMA Spectral Unfolder