ompMC

The original repository is edoerner/ompMC by Edgardo Doerner. This repository is a fork under further development, aimed at integration into the treatment planning toolkits matRad (e0404/matRad) and pyRadPlan (e0404/pyRadPlan).

ompMC is an OpenMP-parallelized, CPU-based Monte Carlo code for coupled photon–electron transport in voxelized geometries. Its physics is a C re-implementation of the EGSnrc condensed-history transport algorithms, restricted to the interactions that matter for megavoltage photon beams, and specialized for the one geometry a treatment planning system needs: a rectilinear dose grid of user-defined materials and densities.

The point of the code is beamlet-based Monte Carlo treatment planning. Rather than producing a single dose distribution, the omc_matrad user code transports histories for many beamlets in one run and returns the dose-influence matrix (Dij) — one sparse column per beamlet — that an optimizer needs for fluence-map optimization. Dose is scored per beamlet and per batch, so a matching variance matrix comes out alongside it. Everything runs on ordinary multi-core CPUs; no GPU, no cluster, no external EGSnrc installation.

Citing ompMC

If you use this code, please cite the work it is based on:

• E. Doerner and P. Caprile, Technical Note: Parallel implementation of the EGSnrc Monte Carlo simulation of ionizing radiation transport using OpenMP, Medical Physics 44(12), 6672–6677 (2017). doi:10.1002/mp.12642

• E. Doerner and P. Caprile, Technical Note: An hybrid parallel implementation for EGSnrc Monte Carlo user codes, Medical Physics 45(8), 3969–3973 (2018). doi:10.1002/mp.13033

• E. Doerner, C. Rebolledo and V. Gomez, Monte Carlo modelling of photon transport using Heterogeneous Computing, Journal of Physics: Conference Series 1043, 012062 (2018). doi:10.1088/1742-6596/1043/1/012062

User codes

Both link against ompmc_core, the transport library built from src/:

• src/ompmc.c — physics: media and PEGS4 data, photon and electron transport, Compton, Rayleigh, pair/triplet, photoelectric, Møller, Bhabha, bremsstrahlung, annihilation, multiple scattering
• src/omc_random.c — random number generation
• src/omc_score.c — dose and variance scoring
• src/omc_utilities.c — input-file parsing and small helpers

Building

See BUILDING.md for the full story — CMake options, how the MATLAB installation is located, and the platform-specific handling of the OpenMP runtime inside a MEX file (which is genuinely fiddly on macOS).

The short version:

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build --config Release --parallel

CMake 3.20+ is required. OpenMP and MATLAB are both optional: without OpenMP the build falls back to serial execution with a warning, and without MATLAB the MEX target is simply skipped. The configure step prints a summary of what was found.

Running omc_dosxyz

-i takes the input file path without the .inp extension; paths inside the input file are resolved relative to the current working directory. From the repository root:

./build/bin/omc_dosxyz -i ucodes/omc_dosxyz/smoke_test -o smoke_test

writes output/smoke_test.3ddose. ucodes/omc_dosxyz/input_file.inp is the fuller example, meant to be run from ucodes/omc_dosxyz/.

Input files are grouped into sections of key = value lines:

# start source definition
mono energy = 20.0
spectrum file = ./../../spectra/mohan6.spectrum
charge = 0
collimator bounds = -2.5 2.5 -2.5 2.5
ssd = 90.0
# end source definition

# start MC control
ncase  = 100000
nbatch = 10
rng seeds = 97 33
# stop MC control

# start geometry
method of input = phantom
phantom file = ./../../phantoms/WATER.egsphant
# stop geometry

# start MC transport
global ecut = 0.521
global pcut = 0.010
pegs file = ./../../pegs4/521icru.pegs4dat
pgs4form file = ./../../pegs4/pgs4form.dat
# stop MC transport

# start VRT
nsplit = 20
# stop VRT

# start ompMC environment
data folder = ./../../data/
output folder = ./../../output/
# stop ompMC environment

Using omc_matrad from MATLAB

addpath('build/bin');
[dij, dijVar] = omc_matrad(cubeRho, cubeMatIx, mcGeo, mcSrc, mcOpt);

• cubeRho — 3D double cube of mass densities
• cubeMatIx — 3D int32 cube of material indices into mcGeo.material
• mcGeo — dose grid: material, xBounds, yBounds, zBounds
• mcSrc — beamlet source: nBixels, iBeam, source position and, per beamlet, the corner and two edge vectors of its aperture
• mcOpt — run settings: nHistories, nBatches, nSplit, charge, global_ecut, global_pcut, randomSeeds, pegsFile, pgs4formFile, dataFolder, outputFolder, and optionally spectrum, spectrumFile, monoEnergy, sourceGeometry ('point' or 'gaussian'), sourceGaussianWidth, relDoseThreshold, verbose, progressCallback, and the variance-reduction keys below

charge picks the source particle: -1 for electrons, 0 for photons, +1 for positrons.

mcOpt.mode — dose-influence matrix or forward dose

'forward_beamlet' computes the dose of a whole weighted field in one go. The collimation is given as one weight per beamlet in mcSrc.bixelWeights, a non-negative vector of length nBixels: a blocked beamlet gets 0, an open one its fluence, a partly transmitting one a fraction of it. This is the fluence map matRad already optimises, so no new geometry is needed.

mcOpt.mode = 'forward_beamlet';
mcSrc.bixelWeights = w;                 % from matRad_fluenceOptimization
[dose, relUnc] = omc_matrad(cubeRho, cubeMatIx, mcGeo, mcSrc, mcOpt);

The result is what dij*w would have been, in Gy for exactly those weights — doubling every weight doubles the dose — but it is reached directly instead of through the matrix. Histories go to the beamlets in proportion to their weight, so a blocked beamlet costs nothing and the run time no longer grows with nBixels.

Two things change meaning in this mode:

• nHistories counts the whole calculation, not one beamlet. Switching a dij run over unchanged therefore divides the statistics by nBixels; multiply it by nBixels to keep them.
• relDoseThreshold does nothing. It prunes columns of a sparse matrix, and there is no matrix here. Note the flip side when comparing the two modes: it is the dij result that is pruned, so set it to 0 for a like-for-like comparison.

relUnc is the relative uncertainty per voxel, 0.9999999 where nothing was deposited — the convention omc_dosxyz writes into a .3ddose file. mcOpt.outputDose = 0 asks for mean deposited energy instead of Gy.

The weights modulate fluence, not spectrum: a leaf transmitting 2% starts 2% of the particles, with the spectrum unhardened. Attenuation in the collimator, its scatter and the beam hardening that goes with it are not modelled. The mode is named for its source model rather than its output, which is what lets 'forward_phsp' sit next to it.

The third output, summary, describes what became of the histories: nHistories, nStarted (how many put a particle into the phantom), nBlocked (how many the collimator stopped) and energyFraction. Mode 'dij' has no equivalent and refuses it — a beamlet that started nothing comes back as a column of zeros, which says so already.

mcOpt.mode = 'forward_phsp' — starting from a phase space

A phase-space file records everything that crossed a plane in an earlier simulation of a treatment head; the sets published at https://www-nds.iaea.org/phsp/ are the output of full models of real linacs. Starting histories from those particles is the difference between modelling the beam and describing it. There is no spectrum and there are no beamlets: the file carries the energy, position and direction of every particle it holds.

mcOpt.mode = 'forward_phsp';
mcSrc.phaseSpace = struct('file', 'Varian_TrueBeam6MV_01');
[dose, relUnc, summary] = omc_matrad(cubeRho, cubeMatIx, mcGeo, mcSrc, mcOpt);

The whole file goes into memory, so it costs about its size on disk — gigabytes for a published dataset. And it was recorded wherever the original simulation scored it, not aimed at your phantom, so it is normal for most histories to start nothing; summary.nStarted is what tells that apart from a transform that is wrong.

One particle per history. A phase space records which particles a single original history left behind, and those are correlated. Drawing them one at a time still gets the dose right on average, but the uncertainty a run reports comes out smaller than the truth by however much they are correlated.

mcSrc.collimator — something in the beam's way

Optional in either forward mode. A transmission grid on a plane, applied by back projection from wherever the source put the particle — so it composes with either source and does not care which side of the plane the particle started on. That is what lets a field be cut out of a phase space recorded above the jaws, as the published ones are.

% A 10 x 10 cm field at the 100 cm isocentre, from a jaw at 40 cm
mcSrc.collimator = struct('z', 40, 'x0', -2, 'y0', -2, 'dx', 4, 'dy', 4, ...
                          'transmission', 1);

With beamlets this is usually unnecessary, the collimation already being in bixelWeights; it earns its place where the weights cannot say what is wanted, such as a block cutting across beamlets or a leaf that transmits.

roulette chooses how a partly transmitting cell is paid for, and both give the same dose in the mean. false multiplies the particle's weight by the fraction and transports it regardless, so a 2% leaf costs a full shower for a fiftieth of the dose — but draws no random numbers at all, leaving every history's random stream where it would have been with the beam open, so a collimated run stays comparable history by history with the open one it came from. true lets the particle through with that probability at full weight instead, spending the time where the dose is, at the price of one random number and more noise per history. Cells that are fully open or fully shut are decided without drawing either way, so an all-or-nothing aperture behaves identically under both.

This is a mask, not a collimator: it attenuates and blocks, but does not scatter and does not harden the spectrum of what it lets through — good for the fluence, poor for the penumbra.

The source spectrum can either be read from a .spectrum file (spectrumFile, default ./spectra/mohan6.spectrum) or passed in directly as mcOpt.spectrum, a struct holding the same information:

mcOpt.spectrum = struct('energy', [1; 2; 3], 'fluence', [0.2; 0.5; 0.3]);

Within a bin the energy is sampled uniformly, as it is for a spectrum read from file.

monoEnergy is the third way: a single kinetic energy in MeV, used for every source particle.

The three are tried in order — spectrum, then spectrumFile, then monoEnergy — and whichever loses is announced rather than silently dropped. Giving none of them uses spectra/mohan6.spectrum.

progressCallback, if given, is a function handle called with a single scalar in [0,1] once per batch and once per finished beamlet; it replaces the built-in waitbar and owns any handle/window lifecycle itself, e.g. mcOpt.progressCallback = @(p) waitbar(p, h, msg);. Without it, a waitbar is shown automatically when verbose >= 2.

Both outputs are sparse, with one column per beamlet and one row per dose-grid voxel; the second output is only computed if requested. Entries below relDoseThreshold (relative to the beamlet maximum) are dropped.

The MEX file calls mexLock() on entry and cannot be unloaded — unloading an OpenMP-using MEX file after a parallel region has run crashes MATLAB. In practice this means a rebuilt MEX file is only picked up after restarting MATLAB. BUILDING.md explains why in detail.

Using ompMC from Python

pip install .

which compiles the extension for the interpreter it is run with; a C++ compiler and a working OpenMP runtime are all it needs. Prebuilt wheels for Linux, macOS and Windows come out of the wheels workflow and carry their own OpenMP runtime, so they need neither.

The wheel bundles the cross section data, PEGS files and spectra, so nothing has to be pointed at the source tree. Four calculations are available, sharing the same phantom and physics:

import numpy as np, ompmc

n = 32
lateral, depth = np.linspace(-8.0, 8.0, n + 1), np.linspace(0.0, 16.0, n + 1)
geometry = ompmc.Geometry(
    lateral, lateral, depth, ["H2O521ICRU"],
    density=np.full((n, n, n), 1.0, order="F"),
    material=np.ones((n, n, n), dtype=np.int32, order="F"),
)

# One dense dose cube from a collimated beam
dose, uncertainty = ompmc.calc_cube(
    geometry,
    ompmc.CollimatedSource(ssd=100.0, x_min=-2, x_max=2, y_min=-2, y_max=2),
    ompmc.Spectrum.monoenergetic(6.0),
    n_histories=100_000, n_batches=10,
)

# ... or one sparse column per beamlet, as scipy.sparse.csc_array
dij = ompmc.calc_dij(geometry, beamlet_source, ompmc.Spectrum.default(),
                     n_histories=100_000, progress=lambda p: print(f"{p:.0%}"))

# ... or the dense cube of a whole weighted field, which is dij @ weights
# computed directly. A blocked beamlet weighs 0 and costs nothing.
dose, uncertainty = ompmc.calc_forward(
    geometry, beamlet_source, weights, ompmc.Spectrum.default(),
    n_histories=100_000,
)

# ... or the same from a linac's own particles, with a 10 x 10 cm field at
# 100 cm cut out of them by a jaw at 40 cm. No spectrum: the file carries one.
dose, uncertainty, summary = ompmc.calc_forward_phsp(
    geometry,
    ompmc.PhaseSpaceSource("Varian_TrueBeam6MV_01"),
    n_histories=1_000_000,
    collimator=ompmc.ApertureMask.rectangle(40.0, -2.0, 2.0, -2.0, 2.0),
)

• Cubes must be Fortran ordered. The transport indexes voxels with the first axis varying fastest, so a C ordered cube would be a silently transposed phantom; it is rejected instead.
• Material indices count from 1, matching matRad's cubeMatIx; 0 means vacuum.
• progress is called with the fraction finished; returning False stops the run, as does Ctrl-C.
• calc_forward is the Python side of mcOpt.mode = 'forward_beamlet' above, with the same two caveats: n_histories counts the whole calculation rather than one beamlet, and the weights modulate fluence rather than spectrum. It takes a collimator= too.
• calc_forward_phsp is the Python side of mcOpt.mode = 'forward_phsp', and carries the same warnings: the whole file goes into memory, one particle starts each history so the reported uncertainty is optimistic, and most histories starting nothing is normal — the file was recorded wherever the original simulation scored it, not aimed at your phantom. It returns a third value, a RunSummary of n_histories, n_started, n_blocked and energy_fraction, which is what tells that apart from a transform that is wrong.
• ompmc.ApertureMask is something in the beam's way, applied by back projection so it composes with either source. roulette=True spends a partly transmitting cell as a survival probability at full weight rather than as a weight multiplier — cheaper behind thick leaves, noisier, and it draws a random number where the default draws none.
• The GIL is released for the whole calculation, so the OpenMP threads run at full speed. The engines keep their state in globals, so one calculation runs at a time per process: use multiprocessing, not threads.
• ompmc.Physics(...) carries the cut-offs, seeds, splitting factor and the variance-reduction keys below.

Variance reduction

Photon transport uses Woodcock (delta) tracking, so photon steps are not stopped at voxel boundaries.

Data files

Tests

Unit tests are built by default (OMPMC_BUILD_TESTS=ON) and registered with CTest:

ctest --test-dir build --output-on-failure

This covers the transport helpers and media data (tests/) plus a short omc_dosxyz smoke run. When the Octave MEX file was built, ctest also drives it through the MEX-side test below. The same test runs unchanged in MATLAB, which needs a MATLAB session:

addpath('build/bin'); addpath('ucodes/omc_matrad');
test_omc_matrad_mex

.github/workflows/build.yml builds and smoke tests every push on Windows x64 (MSVC and MinGW), Linux x64, Linux ARM64, macOS x64 and macOS ARM64.

Changelog

See CHANGELOG.md for a history of changes.

License

GNU General Public License v3.0 — see LICENSE. Copyright (C) 2018-2026 Edgardo Doerner and Niklas Wahl.