Skip to content

Per-file probability normalization constant to enable combining varying statistics simulation output #82

@g5t

Description

@g5t

A possibly-general use case for MCPL output from a simulation is to iteratively build-up a pool of particles which correspond to the same simulation conditions but not necessarily the same overall statistics. For example, @MilanKlausz in #6 seems to have found this necessary for Geant4 simulations, and I have recently been doing the same with McStas simulations.

Details

To summarize the issue:
I wrote an algorithm that repeats a partial-instrument simulation to build-up a list of at-least Ytotal particle states that, e.g., exit a guide for a given chopper-cascade setting, which successively estimates the instrument's transmission, T = Σj Yj / Σj N j, and adjusts the per-iteration initial particle count Nj to reach Σj Yj ≥ Ytotal in the fewest number of iterations, M.

The goal is to then combine the M resulting MCPL files into one file which can be used to feed the next section of the instrument. Since initial particle weights are normalized by the total number of initial particles in the simulation, the particle states from different Nj simulations can not be directly compared or combined.
Instead the weights of all particles recorded in the j-th simulation, must be multiplied by Nj / Σk N k when the M files are combined.

Proposed solutions

All solutions to this general problem likely require the number of initial particles to be recorded alongside the resulting MCPL file.
Since MCPL allows for arbitrary binary 'blob' data, it would be easy to record this number to a specially-named 'blob key'.

Minimal change

If the number of initial particles is recorded to each file, a variant of mctool --merge can perform the per-particle normalization correction when producing the new output file.

I've put together a demonstration of this approach, renormalize, that includes a modified MCPL_output.comp McStas component to add the number of initial particles when writing each file, and a binary tool which mimics mctool --merge and performs the correction.

This solution works, but since this seems to be a general problem an in-MCPL solution is probably better.
To implement the same solution inside of the MCPL distribution, the following decisions/changes should be made:

  • The name of the 'blob key' should be chosen. I've used "mccode_neutron_count", but a more-generic choice likely exists.
  • A method to set and retrieve the number will help MCPL using code. I wrote an encoder and decoder to store the value as a signed 64-bit integer; perhaps a specialized solution for this use case would make more maintainable/useable code?
  • The file-merging function in MCPL should be modified to optionally perform the renormalization and mctool modified to give command line access to the functionality -- perhaps mctool --concat output.mcpl input1.mcpl input2.mcpl ... instead of mctool --merge ....

Fundamental change

If the particle weights recorded in the MCPL file were not normalized by the initial particle count but instead had that normalization applied per-particle while being read back from file, the concatenation of repeated files would not need to modify the (in file) weights.
Instead only the overall normalization would need to be modified.

Effectively, the weights in the 'new' j-th file, $w'$, compared to the 'old' weights, $w$, would be

$$w' = \text{N}_\text{j} w$$

and the combined file would record Σk N k as the number of initial particles.

This may require modifications in all MCPL consuming code, depending on how they write and read particle states.

Calculating the per-particle multiplicative normalization constant (e.g., $1/N$) should be done once as part of opening an mcpl_input_t file. Applying the normalization to each read-particle state is not free, but should hopefully not be a significant hindrance.

Metadata

Metadata

Assignees

No one assigned

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions