pegs4.exe fails silently (exit code 16/18/21, no output) for specific elements/media when using classic MIXT/ELEM card format

[ajilufelix]

I'm generating a custom PEGS4 dataset for ICRP143 paediatric tissue media (57 materials) and running into pegs4.exe failing with no diagnostic output whatsoever — just a nonzero exit code and an empty log.
Setup: Linux64 build, invoked as:
pegs4.exe -i
What works: Standard example inputs (h2o700.pegs4inp, p4imixt.pegs4inp, p4icomp.pegs4inp) all compile fine with exit code 0.
What fails: When I test individual elements as pure ELEM records (classic card format):
ELEM
&INP &END


ENER
&INP AE=0.521,UE=50.511,AP=0.010,UP=50.000 &END
TEST
&INP &END
PWLF
&INP &END
DECK
&INP &END
Results:

C, P → exit 0, valid output (~100-200KB)
H, N, O → exit 21, zero-byte output, zero log output (stdout/stderr both empty)
Mg, Ca → exit 16, zero-byte output, zero log output

No fort.* unit files are created in the working directory on failure either. AE=0.521/UE=50.511/AP=0.010/UP=50.000 is the energy range in all cases.
Questions:

What do exit codes 16/18/21 mean for pegs4.exe?
Why would H/N/O/Mg/Ca fail as pure elements when C/P (and O as part of a mixture in p4imixt/h2o700) succeed?
Is there a required data file or environment setup I might be missing that's element-specific?

[ajilufelix]

Update: I've now ruled out data corruption as the cause. I built a completely independent PEGS4 dataset for the same 57 media using PEGS4GUI (manual entry through the GUI, not the command-line free-format input I mentioned above), and it produces the identical crash — same SIGSEGV, same exact source line, same "monotone = F F" pattern reported by HATCH for all 57 media.
Given two independently-built datasets (one hand-corrected from a buggy source file, one built entirely fresh via PEGS4GUI) both fail identically, this points away from a data-quality issue and toward something specific to this combination:

Source: 100 keV max-energy photon beam (CT localizer-style)
ECUT: 0.521–0.55 MeV (tested both)
Crash: SIGSEGV in electr_, traced to the WHILE ( tuss < range_ep(qel,lelktmp,medium) ) loop in the $COMPUTE-ELOSS-G macro in egsnrc.macros — this loop has no lower bound on lelktmp and appears to walk off the bottom of the range_ep array when a secondary electron's remaining range is smaller than the table's lowest tabulated entry.

Is this a known limitation/edge case with very-low-energy photon sources? Should ECUT be set much higher for this kind of low-energy work, or is there a bounds-check missing in that loop?

[ajilufelix]

Update 2: Confirmed the exact trigger. Setting ECUT=2.0 MeV (far above any physically possible secondary-electron energy for this 100 keV photon source) avoids the crash entirely and the run completes successfully , but at the cost of suppressing electron transport completely (N_step=0, PRESTA-II steps/total = NaN). This confirms the crash is specifically triggered when secondary electrons are actually transported near the low-energy floor of the range table, exactly matching the unguarded WHILE loop I found earlier in $COMPUTE-ELOSS-G.
This isn't a usable workaround for real dosimetry (it suppresses all electron energy deposition), but it does pin down the trigger precisely. I'm going to try intermediate ECUT values to find the lowest setting that avoids the crash while still allowing some electron transport, in case that's useful data too.

[ajilufelix]

Update 4: Ran a full production-scale simulation (10,000,000 histories, ~47 min) at ECUT=0.7 to confirm the smaller test wasn't a fluke. Same result: N_step=0, no electron transport, despite the run completing without crashing. This confirms at scale that ECUT=0.7 avoids the crash only by suppressing electron transport entirely (as expected, since 0.7 MeV total energy / ~189 keV kinetic is well above the ~28 keV Compton-edge maximum for this 100 keV photon source).
I'm now going to test ECUT values closer to the crash boundary (between the working AE floor of 0.521 and the ~0.539 MeV total energy Compton-edge ceiling) to pin down exactly where the SIGSEGV triggers, in case that precision helps diagnose the table-lookup bug.

[ajilufelix]

Tested another proposed workaround — lowering the PEGS4 AE parameter to 0.512 MeV (1 keV kinetic, EGSnrc's absolute minimum) while keeping AE < ECUT (ECUT=0.521 MeV), to give the range table a small buffer below the cutoff rather than having AE and ECUT coincide exactly. Result: this does not crash, but instead hangs indefinitely — 50+ minutes of CPU time at ~100% usage for only 10,000 histories (versus ~47 minutes for a full 10,000,000-history run at a working ECUT), never reaching the first batch print. Confirmed via strace -c -p over a 15-second window: zero system calls observed, indicating a tight, I/O-free computational loop rather than genuine (if slow) progress — consistent with a different known EGSnrc failure mode (#634, "Infinite loop in g application if pegs4 data too low in energy").
So both proposed workarounds tested so far fail, in two different ways:

AE = ECUT = 0.521 → SIGSEGV (reported in Updates 1-3)
AE = 0.512 < ECUT = 0.521 (buffer) → infinite loop, no crash, no output

[rtownson]

I recommend running pegs4 through egs_gui. It has a pegs4 tab, and if you uncheck "ICRU density correction", then you can create pegs4dat files by typing in elements, like you're doing.

Note that pegs4dat files use an approximation of density corrections. A better calculation uses NIST ESTAR, but their website is currently not working properly. The 2026 EGSnrc release will have density corrections from ESTAR directly integrated, so you will be able to define materials right in your input file and get the most accurate results.

H, N, and O need to be diatomic molecules, so pegs4 is meant to fail with an error message telling you this. It works for me when I run it through egs_gui. Mg and Ca should work, so let me know if you also have issues with them when using egs_gui. If it also fails for you, then it might be an issue with newer compilers.

[ajilufelix]

Peg4 has been run through gui successfully but when I set ECUT = 0.521 it crashed and 512 never crashed but stucked At ECUT = 0.7 MeV — Completed, but No Electron Transport
The simulation completed successfully through all batches. However, post-run diagnostics reported N_step = 0 — zero electron transport steps occurred anywhere in the run.

This is a direct physical consequence of the source energy: the maximum possible secondary-electron kinetic energy from Compton scattering of a 100 keV photon (the Compton edge) is approximately 28 keV — well below the ≈189 keV kinetic-energy threshold implied by ECUT = 0.7 MeV total energy. Every secondary electron therefore had its energy deposited locally at the point of creation rather than being transported through the medium.

Consequence: the resulting dose distribution reflects photon interaction dose only. It does not model the spatial spread of electron energy deposition, which can be significant near tissue and organ boundaries.

[ajilufelix]

Subject: DOSXYZnrc — correct global energy balance but per-voxel .3ddose values off by ~20 orders of magnitude (phase-space source, isource=2)

Hi all,

I'm running into a persistent issue with DOSXYZnrc that I can't resolve
through input-file configuration, and I'd appreciate any insight.

Setup

• DOSXYZnrc simulation on a CT-derived voxel phantom (ICRP-143 reference
  phantom family, converted to .egsphant)
• Source: isource=2 (full phase-space source), generated by a custom
  BEAMnrc accelerator model (kV diagnostic imaging tube: tungsten target,
  Ag/Cu filter, bowtie filter, collimator)
• Phase-space file independently verified clean via BEAMDP: 100% photons,
  sensible energy range (0.010–0.119 MeV), no NaN/degenerate records
• PEGS4 dataset: ICRP-143 tissue media, AE=0.521 MeV (appropriate for the
  ~100 keV beam energy)
• NCASE = 10,000 (phase-space source, with recycling as needed)

Symptom

DOSXYZnrc's own global energy-balance report is physically sensible:

Fraction of incident energy deposited in the phantom = 0.3694 (36.9%)

But the resulting .3ddose file (and the raw per-voxel depth-dose tables
printed to .egslst with doseprint=1) show almost no real per-voxel dose:

Total voxels:              55,123,640
Nonzero-dose voxels:       1
Max dose:                  6.5e-34 Gy/particle
Every printed value carries a reported uncertainty of exactly -99.9%

For comparison, a validated water-phantom test case (much smaller grid,
~10,000 voxels, different source) run earlier in this project produced
sensible per-voxel dose values in the expected 1e-11 to 1e-14 Gy/particle
range with normal (<50%) uncertainties — so I know a "healthy" run should
look nothing like this.

What I've ruled out so far

1. Custom phantom-cropping — ran the identical setup against the original,
   uncropped full-body phantom; identical anomaly present.
2. doseprint=0 suppressing output — set doseprint=1, confirmed genuine
   detailed per-voxel output is being printed (17M-line .egslst log); the
   values themselves, not their suppression, are the problem.
3. Beam mis-centered on the phantom's bounding-box corner rather than its
   anatomical center — found and fixed this real bug (confirmed via a 12x
   increase in global energy-deposition fraction, 2.97% -> 36.26%), but it
   did not resolve the per-voxel magnitude anomaly.
4. Phase-space source file quality/composition — independently verified
   clean via BEAMDP.
5. Dose normalization (NP, per the isource=2 formula in Section 13 of the
   Users Manual) — manually recomputed from NINCSRC/nshist/NCASE, matches
   DOSXYZnrc's own printed value almost exactly (605,591 vs my calculated
   605,495). This is a real, legitimate ~60x normalization factor, but
   cannot explain a ~1e20 magnitude gap.
6. phicol (source plane rotation angle) left at the default 0 instead of
   the documented required value of 180 for BEAMnrc-generated phase-space
   sources (Section 4.4) — corrected to 180, confirmed via the run's own
   echoed "Source plane rotation angle: 180.0000 degrees" — did not change
   the anomaly at all.
7. Confirmed this anomaly was present in the very first DOSXYZnrc run at
   the start of this project (different phantom, different beam
   coordinates), so it predates every change made while debugging it.

Question

Given correct global energy accounting but per-voxel dose values off by
~20 orders of magnitude and universally flagged with -99.9% uncertainty,
does this match any known issue — e.g. a specific compiler/build problem,
a known bug in a particular DOSXYZnrc/EGSnrc version, or a documented
gotcha with isource=2 phase-space sources that isn't covered in Section
4.4 of the Users Manual?

I'm currently planning to break with gdb directly at the dose-accumulation
line in the transport loop to see whether the deposited-energy value
itself is already wrong at the point of accumulation, or whether it's
corrupted only during the final array write — but if this is a known
issue, I'd very much like to skip that step.

System details:
EGSnrc version 2025a
Ubuntu version 22.04
gfortran 11.4.0

Thanks in advance,
Felix

[ajilufelix]

Subject: BEAMnrc/DOSXYZnrc — phase-space particle weights show geometric decay in powers of (1/n_split), consistent with recycling compounding the splitting weight-reduction

Hi all,

I've tracked down a specific, reproducible issue in a BEAMnrc-generated
phase-space file that I believe explains a severe downstream dose-magnitude
anomaly in DOSXYZnrc, and I'd like to check whether this is a known issue
or a misconfiguration on my end.

Setup

• Custom BEAMnrc accelerator (kV diagnostic imaging tube: tungsten target,
  Ag/Cu filter, bowtie filter, collimator), 100 kVp
• Bremsstrahlung splitting enabled: NBRSPL = 500 (uniform)
• Phase-space file scored at the collimator output (CM4), MODE0 format
• This file subsequently used as a DOSXYZnrc isource=2 source, with
  recycling (NRCYCL) as needed to reach the requested NCASE

Symptom (downstream, in DOSXYZnrc)

DOSXYZnrc's global energy-balance accounting is physically sensible
(~36% of incident energy deposited in the phantom), but the resulting
per-voxel .3ddose values are essentially all zero -- with only a
statistically negligible fraction of voxels showing any dose at all, and
those at implausibly small magnitudes (down to ~1e-34 Gy/particle, versus
an expected ~1e-11 to 1e-14 Gy/particle for a validated reference case).

Root cause, isolated via gdb + direct binary inspection

Breaking directly at the phase-space read statement in DOSXYZnrc
(dosxyznrc_linux64.F, around the weight=WT_PHSP_SHORT assignment)
showed wildly inconsistent particle weights between consecutive reads
(e.g. 6.4e-18, then 1.9e-13, then values differing by many orders of
magnitude with no apparent pattern at first).

I then wrote a script to scan the raw phase-space file directly
(bypassing DOSXYZnrc entirely) and compute the full weight distribution
across all ~165,000 particle records. Key statistics:

Max weight:    2.000000e-03   (exactly 1/500 -- matches NBRSPL=500)
Min weight:    2.220446e-19
Mean weight:   3.461321e-04
Median weight: 2.155056e-05
39.1% of records have |weight| < 1e-6

Critically, the minimum observed weight (2.22e-19) is very close to
(1/500)^7 = 1.28e-19, and the median (2.16e-05) is close to
(1/500)^2 = 4.00e-06. The weight distribution appears to follow a
geometric progression in powers of (1/n_split), rather than clustering
near a single expected value of 1/500 for all split particles.

My interpretation

This pattern is consistent with the 1/n_split weight reduction being
applied repeatedly to the same particle across multiple splits/recycles
during phase-space generation or subsequent recycling, rather than being
applied once per original interaction and then preserved across reuse.
Only the never-recycled fraction of particles retains a "normal" weight
near 1/500; particles reused N times appear to carry weight
(1/500)^N, contributing almost nothing to any dose tally, which explains
both the near-zero per-voxel .3ddose values and the fact that global
energy-balance accounting (which is presumably computed independently of
per-voxel weighting) still looks correct.

Question

Does this match a known issue with combining bremsstrahlung splitting and
phase-space recycling (either within a single BEAMnrc run, or across
DOSXYZnrc's NRCYCL-based reuse of a phase-space source)? Is there a
documented requirement to disable one when using the other, or a specific
input combination known to cause this kind of weight-compounding?

I'm happy to share the phase-space file, the full .egsinp for both the
BEAMnrc accelerator and the DOSXYZnrc simulation, and the small Python
script used to compute the weight distribution above, if useful.

System details:

• EGSnrc git commit: 9edee3e (2025-12-04 13:26:33 -0500)
• OS: Ubuntu 22.04.5 LTS
• Compiler: GNU Fortran (Ubuntu 11.4.0-1ubuntu1~22.04.3) 11.4.0

Thanks in advance,
Felix

---

Appendix: ruled out before reaching this diagnosis (in case relevant):

• Custom phantom-cropping tool -- confirmed not the cause (identical
  anomaly present with the original, uncropped full-body phantom)
• doseprint=0 suppressing output -- ruled out (doseprint=1 confirmed
  genuine detailed per-voxel output is generated; magnitude anomaly is
  present in that raw output too, not just the binary .3ddose file)
• Beam mis-centered on the phantom's bounding-box corner rather than its
  anatomical center -- real, separate bug, found and fixed (confirmed via
  a 12x increase in global energy-deposition fraction), but did not
  resolve this weight/dose-magnitude issue
• phicol left at the default 0 instead of the documented required value
  of 180 for BEAMnrc-generated phase-space sources (Users Manual Section
  4.4) -- corrected, confirmed via echoed "Source plane rotation angle:
  180.0000 degrees", did not change the anomaly
• Dose normalization (NP, per the isource=2 formula in Users Manual
  Section 13) -- manually recomputed from NINCSRC/nshist/NCASE, matches
  DOSXYZnrc's own printed value almost exactly
• PHSP_RECL (direct-access record length used to open the phase-space
  file) -- confirmed correct (28 bytes, matching the file's real MODE0
  format) via gdb
• ismode (MODE0 vs MODE2 record format flag) -- confirmed correctly set
  to 0, matching the file's real format

[rtownson]

Make sure 'howfarless' is off - it is only valid for entirely homogenous phantoms.

By default, DOSXYZnrc zeros doses with >50% uncertainty. Since you're only running 10k particles, this is likely the issue. This was originally implemented to stop high uncertainty doses from throwing off the reporting of the "maximum dose".

In $EGS_HOME/dosxyznrc/dosxyznrc_user_macros.mortran, change line 105 from:

REPLACE {$DOSEZERO} WITH {1};

to instead be:

REPLACE {$DOSEZERO} WITH {0};

Then recompile dosxyznrc. You can do this using the dosxyz_gui, or on the command line from within the egs_home/dosxyznrc folder with make clean; make on linux, or mingw32-make clean and then mingw32-make on windows.

[ajilufelix]

Thank you for the reply

[ajilufelix]

Subject: BEAMnrc — phase-space particle weights show geometric decay in powers of (1/n_split), specifically tied to multi-batch runs (single-batch runs are clean)

Hi all,

Follow-up to my earlier post (below) with a much more precise finding after
further testing.

Refined diagnosis

I initially suspected combining bremsstrahlung splitting (NBRSPL=500) with
photon forcing was the cause. I turned photon forcing off and reran with a
small NCASE (100,000 histories) that completed in a single batch --
this produced a genuinely clean phase-space file: every particle had
identical weight, exactly 1/500 (0.002), confirmed via both BEAMnrc's own
run-summary (MIN=MAX=0.0020000) and an independent script that scans every
raw record in the file.

However, rerunning the same accelerator (photon forcing still off,
NBRSPL=500 unchanged) with a larger NCASE=10,000,000 -- which necessarily
spans 10 separate batches -- reproduced the exact same corrupted weight
distribution as before:

Records read: 165,128 (header claims only 28,549)
Min weight:    2.220446e-19
Max weight:    2.000000e-03   (correct, = 1/500)
32.3% of records have |weight| < 1e-6

The minimum observed weight (2.22e-19) is close to (1/500)^7 = 1.28e-19.
Critically, with 10 batches total, if every particle were compounding
once per batch we'd expect values as small as (1/500)^10 -- several more
orders of magnitude smaller than observed. Instead it looks like only a
subset of particles are affected, each compounded roughly 7 times, while
the majority retain the correct, uniform 1/500 weight.

This points to something specific to carrying particles across batch
boundaries during a multi-batch run (rather than forcing+splitting in
general) -- e.g. some particles from one batch's phase-space-writing state
being incorrectly re-processed (and re-split) in a subsequent batch.

Question

Does this match a known issue with phase-space output across multiple
NCASE batches when combined with bremsstrahlung splitting? Is there a
batch-related state (e.g. related to IPARALLEL, or internal
phase-space-writing counters) that isn't being correctly reset/tracked
between batches in this scenario?

As a workaround, restricting a run to complete within a single batch (by
keeping NCASE low enough) produces clean, correctly-weighted output --
but this obviously limits the statistics achievable in one run.

System details:

• EGSnrc git commit: 9edee3e (2025-12-04 13:26:33 -0500)
• OS: Ubuntu 22.04.5 LTS
• Compiler: GNU Fortran (Ubuntu 11.4.0-1ubuntu1~22.04.3) 11.4.0

Thanks again for the help so far -- the ihowfarless and $DOSEZERO pointers
were both spot-on and much appreciated.

Felix

---

Original post:

Subject: BEAMnrc/DOSXYZnrc — phase-space particle weights show geometric decay in powers of (1/n_split), consistent with recycling compounding the splitting weight-reduction

Hi all,

I've tracked down a specific, reproducible issue in a BEAMnrc-generated
phase-space file that I believe explains a severe downstream dose-magnitude
anomaly in DOSXYZnrc, and I'd like to check whether this is a known issue
or a misconfiguration on my end.

Setup

• Custom BEAMnrc accelerator (kV diagnostic imaging tube: tungsten target,
  Ag/Cu filter, bowtie filter, collimator), 100 kVp
• Bremsstrahlung splitting enabled: NBRSPL = 500 (uniform)
• Phase-space file scored at the collimator output (CM4), MODE0 format
• This file subsequently used as a DOSXYZnrc isource=2 source, with
  recycling (NRCYCL) as needed to reach the requested NCASE

Symptom (downstream, in DOSXYZnrc)

DOSXYZnrc's global energy-balance accounting is physically sensible
(~36% of incident energy deposited in the phantom), but the resulting
per-voxel .3ddose values are essentially all zero -- with only a
statistically negligible fraction of voxels showing any dose at all, and
those at implausibly small magnitudes (down to ~1e-34 Gy/particle, versus
an expected ~1e-11 to 1e-14 Gy/particle for a validated reference case).

Root cause, isolated via gdb + direct binary inspection

Breaking directly at the phase-space read statement in DOSXYZnrc
(dosxyznrc_linux64.F, around the weight=WT_PHSP_SHORT assignment)
showed wildly inconsistent particle weights between consecutive reads
(e.g. 6.4e-18, then 1.9e-13, then values differing by many orders of
magnitude with no apparent pattern at first).

I then wrote a script to scan the raw phase-space file directly
(bypassing DOSXYZnrc entirely) and compute the full weight distribution
across all ~165,000 particle records. Key statistics:

Max weight:    2.000000e-03   (exactly 1/500 -- matches NBRSPL=500)
Min weight:    2.220446e-19
Mean weight:   3.461321e-04
Median weight: 2.155056e-05
39.1% of records have |weight| < 1e-6

Critically, the minimum observed weight (2.22e-19) is very close to
(1/500)^7 = 1.28e-19, and the median (2.16e-05) is close to
(1/500)^2 = 4.00e-06. The weight distribution appears to follow a
geometric progression in powers of (1/n_split), rather than clustering
near a single expected value of 1/500 for all split particles.

My interpretation

This pattern is consistent with the 1/n_split weight reduction being
applied repeatedly to the same particle across multiple splits/recycles
during phase-space generation or subsequent recycling, rather than being
applied once per original interaction and then preserved across reuse.
Only the never-recycled fraction of particles retains a "normal" weight
near 1/500; particles reused N times appear to carry weight
(1/500)^N, contributing almost nothing to any dose tally, which explains
both the near-zero per-voxel .3ddose values and the fact that global
energy-balance accounting (which is presumably computed independently of
per-voxel weighting) still looks correct.

Question

Does this match a known issue with combining bremsstrahlung splitting and
phase-space recycling (either within a single BEAMnrc run, or across
DOSXYZnrc's NRCYCL-based reuse of a phase-space source)? Is there a
documented requirement to disable one when using the other, or a specific
input combination known to cause this kind of weight-compounding?

I'm happy to share the phase-space file, the full .egsinp for both the
BEAMnrc accelerator and the DOSXYZnrc simulation, and the small Python
script used to compute the weight distribution above, if useful.

System details:

• EGSnrc git commit: 9edee3e (2025-12-04 13:26:33 -0500)
• OS: Ubuntu 22.04.5 LTS
• Compiler: GNU Fortran (Ubuntu 11.4.0-1ubuntu1~22.04.3) 11.4.0

Thanks in advance
Felix

---

Appendix: ruled out before reaching this diagnosis (in case relevant):

• Custom phantom-cropping tool -- confirmed not the cause (identical
  anomaly present with the original, uncropped full-body phantom)
• doseprint=0 suppressing output -- ruled out (doseprint=1 confirmed
  genuine detailed per-voxel output is generated; magnitude anomaly is
  present in that raw output too, not just the binary .3ddose file)
• Beam mis-centered on the phantom's bounding-box corner rather than its
  anatomical center -- real, separate bug, found and fixed (confirmed via
  a 12x increase in global energy-deposition fraction), but did not
  resolve this weight/dose-magnitude issue
• phicol left at the default 0 instead of the documented required value
  of 180 for BEAMnrc-generated phase-space sources (Users Manual Section
  4.4) -- corrected, confirmed via echoed "Source plane rotation angle:
  180.0000 degrees", did not change the anomaly
• Dose normalization (NP, per the isource=2 formula in Users Manual
  Section 13) -- manually recomputed from NINCSRC/nshist/NCASE, matches
  DOSXYZnrc's own printed value almost exactly
• PHSP_RECL (direct-access record length used to open the phase-space
  file) -- confirmed correct (28 bytes, matching the file's real MODE0
  format) via gdb
• ismode (MODE0 vs MODE2 record format flag) -- confirmed correctly set
  to 0, matching the file's real format

[ajilufelix]

FINAL UPDATE — Resolved

False alarm on both the "photon forcing" and "batch boundary" theories above
-- the actual cause was much simpler and entirely on my end: the
phase-space output file was never being deleted/truncated between runs.
Each new run's file only had its header metadata updated to reflect the
latest run's particle count, while the underlying binary records were
being silently appended to (or otherwise not cleared), accumulating
corrupted leftover data from an much earlier test run (my very first
attempt with photon forcing "on", which itself may have produced bad
output for unrelated reasons I never fully diagnosed).

Explicitly deleting the .egsphsp1 file before each fresh run immediately
resolved everything: file size, record count, and weight distribution all
now line up exactly as expected (uniform weight = 1/NBRSPL for every
particle, header particle count matching the actual record count exactly).

Apologies for the noise on the two intermediate theories -- leaving them
in below in case anyone finds the diagnostic approach (gdb + direct binary
record inspection + full weight-distribution scanning) useful for a
similar-looking problem that turns out to have a different actual cause.

Thanks again for the ihowfarless and $DOSEZERO help, both of which remain
valid, useful fixes independent of this file-handling issue.

---

Update 2:

Subject: BEAMnrc — phase-space particle weights show geometric decay in powers of (1/n_split), specifically tied to multi-batch runs (single-batch runs are clean)

Hi all,

Follow-up to my earlier post (below) with a much more precise finding after
further testing.

Refined diagnosis

I initially suspected combining bremsstrahlung splitting (NBRSPL=500) with
photon forcing was the cause. I turned photon forcing off and reran with a
small NCASE (100,000 histories) that completed in a single batch --
this produced a genuinely clean phase-space file: every particle had
identical weight, exactly 1/500 (0.002), confirmed via both BEAMnrc's own
run-summary (MIN=MAX=0.0020000) and an independent script that scans every
raw record in the file.

However, rerunning the same accelerator (photon forcing still off,
NBRSPL=500 unchanged) with a larger NCASE=10,000,000 -- which necessarily
spans 10 separate batches -- reproduced the exact same corrupted weight
distribution as before:

Records read: 165,128 (header claims only 28,549)
Min weight:    2.220446e-19
Max weight:    2.000000e-03   (correct, = 1/500)
32.3% of records have |weight| < 1e-6

The minimum observed weight (2.22e-19) is close to (1/500)^7 = 1.28e-19.
Critically, with 10 batches total, if every particle were compounding
once per batch we'd expect values as small as (1/500)^10 -- several more
orders of magnitude smaller than observed. Instead it looks like only a
subset of particles are affected, each compounded roughly 7 times, while
the majority retain the correct, uniform 1/500 weight.

This points to something specific to carrying particles across batch
boundaries during a multi-batch run (rather than forcing+splitting in
general) -- e.g. some particles from one batch's phase-space-writing state
being incorrectly re-processed (and re-split) in a subsequent batch.

Question

Does this match a known issue with phase-space output across multiple
NCASE batches when combined with bremsstrahlung splitting? Is there a
batch-related state (e.g. related to IPARALLEL, or internal
phase-space-writing counters) that isn't being correctly reset/tracked
between batches in this scenario?

As a workaround, restricting a run to complete within a single batch (by
keeping NCASE low enough) produces clean, correctly-weighted output --
but this obviously limits the statistics achievable in one run.

System details:

• EGSnrc git commit: 9edee3e (2025-12-04 13:26:33 -0500)
• OS: Ubuntu 22.04.5 LTS
• Compiler: GNU Fortran (Ubuntu 11.4.0-1ubuntu1~22.04.3) 11.4.0

Thanks again for the help so far -- the ihowfarless and $DOSEZERO pointers
were both spot-on and much appreciated.

Felix

---

Original post:

Subject: BEAMnrc/DOSXYZnrc — phase-space particle weights show geometric decay in powers of (1/n_split), consistent with recycling compounding the splitting weight-reduction

Hi all,

I've tracked down a specific, reproducible issue in a BEAMnrc-generated
phase-space file that I believe explains a severe downstream dose-magnitude
anomaly in DOSXYZnrc, and I'd like to check whether this is a known issue
or a misconfiguration on my end.

Setup

• Custom BEAMnrc accelerator (kV diagnostic imaging tube: tungsten target,
  Ag/Cu filter, bowtie filter, collimator), 100 kVp
• Bremsstrahlung splitting enabled: NBRSPL = 500 (uniform)
• Phase-space file scored at the collimator output (CM4), MODE0 format
• This file subsequently used as a DOSXYZnrc isource=2 source, with
  recycling (NRCYCL) as needed to reach the requested NCASE

Symptom (downstream, in DOSXYZnrc)

DOSXYZnrc's global energy-balance accounting is physically sensible
(~36% of incident energy deposited in the phantom), but the resulting
per-voxel .3ddose values are essentially all zero -- with only a
statistically negligible fraction of voxels showing any dose at all, and
those at implausibly small magnitudes (down to ~1e-34 Gy/particle, versus
an expected ~1e-11 to 1e-14 Gy/particle for a validated reference case).

Root cause, isolated via gdb + direct binary inspection

Breaking directly at the phase-space read statement in DOSXYZnrc
(dosxyznrc_linux64.F, around the weight=WT_PHSP_SHORT assignment)
showed wildly inconsistent particle weights between consecutive reads
(e.g. 6.4e-18, then 1.9e-13, then values differing by many orders of
magnitude with no apparent pattern at first).

I then wrote a script to scan the raw phase-space file directly
(bypassing DOSXYZnrc entirely) and compute the full weight distribution
across all ~165,000 particle records. Key statistics:

Max weight:    2.000000e-03   (exactly 1/500 -- matches NBRSPL=500)
Min weight:    2.220446e-19
Mean weight:   3.461321e-04
Median weight: 2.155056e-05
39.1% of records have |weight| < 1e-6

Critically, the minimum observed weight (2.22e-19) is very close to
(1/500)^7 = 1.28e-19, and the median (2.16e-05) is close to
(1/500)^2 = 4.00e-06. The weight distribution appears to follow a
geometric progression in powers of (1/n_split), rather than clustering
near a single expected value of 1/500 for all split particles.

My interpretation

This pattern is consistent with the 1/n_split weight reduction being
applied repeatedly to the same particle across multiple splits/recycles
during phase-space generation or subsequent recycling, rather than being
applied once per original interaction and then preserved across reuse.
Only the never-recycled fraction of particles retains a "normal" weight
near 1/500; particles reused N times appear to carry weight
(1/500)^N, contributing almost nothing to any dose tally, which explains
both the near-zero per-voxel .3ddose values and the fact that global
energy-balance accounting (which is presumably computed independently of
per-voxel weighting) still looks correct.

Question

Does this match a known issue with combining bremsstrahlung splitting and
phase-space recycling (either within a single BEAMnrc run, or across
DOSXYZnrc's NRCYCL-based reuse of a phase-space source)? Is there a
documented requirement to disable one when using the other, or a specific
input combination known to cause this kind of weight-compounding?

I'm happy to share the phase-space file, the full .egsinp for both the
BEAMnrc accelerator and the DOSXYZnrc simulation, and the small Python
script used to compute the weight distribution above, if useful.

System details:

• EGSnrc git commit: 9edee3e (2025-12-04 13:26:33 -0500)
• OS: Ubuntu 22.04.5 LTS
• Compiler: GNU Fortran (Ubuntu 11.4.0-1ubuntu1~22.04.3) 11.4.0

Thanks in advance,
Felix

---

Appendix: ruled out before reaching this diagnosis (in case relevant):

• Custom phantom-cropping tool -- confirmed not the cause (identical
  anomaly present with the original, uncropped full-body phantom)
• doseprint=0 suppressing output -- ruled out (doseprint=1 confirmed
  genuine detailed per-voxel output is generated; magnitude anomaly is
  present in that raw output too, not just the binary .3ddose file)
• Beam mis-centered on the phantom's bounding-box corner rather than its
  anatomical center -- real, separate bug, found and fixed (confirmed via
  a 12x increase in global energy-deposition fraction), but did not
  resolve this weight/dose-magnitude issue
• phicol left at the default 0 instead of the documented required value
  of 180 for BEAMnrc-generated phase-space sources (Users Manual Section
  4.4) -- corrected, confirmed via echoed "Source plane rotation angle:
  180.0000 degrees", did not change the anomaly
• Dose normalization (NP, per the isource=2 formula in Users Manual
  Section 13) -- manually recomputed from NINCSRC/nshist/NCASE, matches
  DOSXYZnrc's own printed value almost exactly
• PHSP_RECL (direct-access record length used to open the phase-space
  file) -- confirmed correct (28 bytes, matching the file's real MODE0
  format) via gdb
• ismode (MODE0 vs MODE2 record format flag) -- confirmed correctly set
  to 0, matching the file's real format