Repository navigation
MCPLGen: normalise the weights by stat:sum:initial_ray_count when the file has it (McStas MCPL_output weight_mode=1) #61
Description
Claude told me that we will need special handling of the MCPL files that are comming from the latest McStas versions due to a change in what weights are stored by a new default in weight_mode. It think this new default breaks any Geant4 simulation that uses MCPLGen, so it should be handled within the framework rather than in the user codes.
The following is an AI written summary of the issue and high level suggestions about how it could be handled.
I suggest adding an MCPLGen option that divides the particle weights by stat:sum:initial_ray_count when the file provides it, as files written with the new McStas MCPL_output default (weight_mode=1) do. It should be on by default, so that the results of existing code stay the same with new McStas versions.
Background. Since McStas 3.10.0 (McCode commit dacc4b87), MCPL_output writes by default with weight_mode=1: each particle weight is stored multiplied by the number of rays of the run (ncount), and the file header gets stat:sum:initial_ray_count = ncount. The physical weight (n/s for a source in n/s) is the stored weight divided by initial_ray_count. Earlier McStas versions (and weight_mode=2) stored the physical weight directly, without initial_ray_count. The weight_mode parameter exists since May 2025 (MCPL >= 2.1); only the default changed. weight_mode=2 now prints a warning that it "will soon become unavailable" (McCode discussion #2076), so files with initial_ray_count will become the only kind McStas writes.
What MCPLGen does now. It passes the stored weight to the primary vertex as it is (mod.cc#L272-L281), and the stat:sum values are not used anywhere. So, with files from a current McStas, the weights in Geant4 (and in every analysis output: histograms, GRIFF, MCPL written by MCPLWriter) are ncount times larger than with the same simulation from an older McStas, and two McStas runs with different ray counts give different results. Nothing on the reading side (MCPL C or Python API) rescales the weights either.
Proposal. A boolean MCPLGen parameter, named after the MCPL header value it uses:
//Divide the weights by stat:sum:initial_ray_count if the file has it (e.g. McStas MCPL_output, weight_mode=1),
//so that the weights are per initial ray of the simulation that wrote the file and do not depend on its ray count:
addParameterBoolean("normalise_by_initial_ray_count",true);
How I think the parameter should be handled:
- When the file is opened, MCPLGen reads the value once with
mcpl_hdr_stat_sum(f,"initial_ray_count")(the particles are not read for this) and prints whether it divides the weights by it. - Parameter true (the default):
- value present and positive: every weight is divided by it, also weights of 1 and a universal weight (the current code skips
SetWeightfor weights equal to 1, so the division has to come before that check); - value not present (-2: older McStas files,
weight_mode=2, files from other codes): the weights are used as stored, exactly as now; - value present but not finalised (-1, e.g. an interrupted McStas run), or zero or negative: stop with an error that names the parameter, since such a file cannot be normalised; setting the parameter to false reads it with the stored weights.
- value present and positive: every weight is divided by it, also weights of 1 and a universal weight (the current code skips
- Parameter false: the weights as stored, for every file (the current behaviour).
- Partial reading (
skip_events, a limited number of events): the division always uses the value of the whole file, so the jobs that each process a part of a file stay additive (the sum of their outputs is the result of the whole file). - Merged files:
mcpltool --mergeadds up theinitial_ray_countof the merged files, so repeated McStas runs merged before Geant4 normalise correctly. Files that are not repeats of the same simulation (e.g. different McStas outputs of one run) are read separately, each with its own value. - Output files: MCPL files written by MCPLWriter carry the normalised weights and do not copy the
initial_ray_countof the input (MCPLWriter does not transfer stat:sum values now, and should not start to), so a file written by Geant4 and read again by MCPLGen is not divided a second time.
Activity
tkittel commentedon Nov 5, 2026
Right...
Question: Should dgcode's MCPLWriter also start to record the number of geant4 events in a separate stat:sum? Should it somehow transfer initial_ray_count?
tkittel commentedon Nov 6, 2026
@MilanKlausz have a look at #62
MilanKlausz commentedon Nov 6, 2026
Thanks, I tried the PR on a real McStas → Geant4 chain and the normalisation itself works fine. Claude run into one real bug that probably only affect Intel MACs, and a few smaller things. See the AI summary below:
Tested on 9504d6d with real McStas 3.10.1 files: results match
normalise_by_initial_ray_count=0divided by ncount, and McStas weights of 1 n/s give Geant4 weights of exactly 1. To your question: yes togeant4_nevts, and transferringinitial_ray_countworks for whole-file jobs (-jand merged repeats).1. Bug on Intel macOS: the transfer silently fails and
testmcplstatsumsfails. Thedynamic_casttoInitialRayCountInfoinMCPLWriterreturns null, so the output has noinitial_ray_count. The class is header-only, so each pymod has its own typeinfo; Python imports withRTLD_LOCAL, so the copies aren't merged, and libc++ on x86_64 Apple compares typeinfo by address. Linux and arm64 macOS compare names, so CI doesn't see it.Tested fix: make the destructor non-inline and define it in
G4MCPL/libsrc.2. Jobs split with
skip_eventslose the normalisation when force-merged (tested with fix 1). Each output keeps the weights multiplied by ncount, with the value unavailable.mcpltool --mergerefuses the parts because the generator settings differ.--forcemergeworks but drops all stat:sum values, so reading the result gives no error and weights too large by ncount (initial_geant4_weights= 5000 instead of 5). I have no good fix: keeping Geant4 weights forskip_events != 0fails because the first job hasskip_events=0. A note in the docs may be enough.3.
initial_geant4_weightsdepends on the flag, but the file doesn't. Output weights are 1000 either way, while the value is 5 or 5000. A sentence in the docs would do.4. Correction to the issue and the MCPLGen error message: an interrupted McStas run only gives -1 on a hard kill. On SIGTERM or Ctrl-C, McStas 3.10.1 records the requested ray count, so the file looks complete and the intensity comes out too low (0.073 instead of 628 n/s in a test). This should be raised on the McStas side and needs no change in the PR.
tkittel commentedon Nov 7, 2026
@MilanKlausz So I just merged #62 with what I think is a suitable fix for this. It won't be auto-magic for all use-cases, since that is not really possible (G4 does not have 1-particle-per-event always).
Some details from Claude:
Intel macOS: MCPLWriter no longer uses a dynamic_cast to find the information from MCPLGen. Instead there is now a small global registry in the G4MCPL library (G4MCPL/OutputStatSums.hh), which MCPLGen uses to provide the stat:sum:initial_ray_count of its input file. The ExprParser exception classes, which had the same problem when caught in Python modules, now have out-of-line destructors. CI now also runs on Intel macOS (macos-15-intel), and that job passes. A retest on your Intel Mac would still be very welcome.
Split jobs: we concluded that the initial_ray_count convention assumes that each particle is an independent ray, which is usually not true for files written by Geant4 (a single event can give several correlated particles). So by default, MCPLWriter no longer claims the convention: the output files have the weights as used in Geant4, and the input file's value is only recorded as stat:sum:original_initial_ray_count, for reference. Outputs of split jobs then simply add up, also after --forcemerge. Users who want the convention in the output can request it explicitly with the opt_initial_ray_count keyword of --mcpl (or setWriteInitialRayCount(True)): the value is then recorded as stat:sum:initial_ray_count, with the weights in the same convention as in the input file.
stat:sum:initial_geant4_weights is now always in the same convention as the weights in the file.
Interrupted McStas runs: reported on the McCode side.
When multiprocessing, values describing the whole job are now divided between the processes, so the merged file gets the value. Setting them only in the merged file would need MCPL to support updating stat:sum values of existing files (mctools/mcpl#254, followed up in #63).
assisted-by: claude
tkittel commentedon Nov 9, 2026
Update on where everything stands now. In short: everything that can be addressed before the next mcpl release has been, and new simplebuild and dgcode releases are out already.
Claude says:
Update: everything from this issue is now released in simple-build-dgcode 2.1.0, on PyPI and conda-forge.
Since my last comment:
stat:sum:initial_ray_count = -1now only says what -1 means in MCPL (the value is not available), without guessing at the cause.sb_mcplextracommands, including invalid ones. The error handling no longer depends on type information being shared between libraries, which is what failed on Intel macOS. CI now also runs on Intel macOS.Follow-ups elsewhere:
tkittel commentedon Nov 9, 2026
@MilanKlausz is that good enough to close here, or is anything still blocking you in the short term?
MilanKlausz commentedon Nov 9, 2026
@tkittel yes, this covers everythin I need now, thank you!
One small thing about your comment from 7 Nov (that Claude noticed):
"Users who want the convention in the output can request it explicitly with the opt_initial_ray_count keyword of --mcpl (or setWriteInitialRayCount(True))"
It calls the new keyword of --mcpl
opt_initial_ray_count, but the code and the docs call itopt_irc.