MCPL >= 2.2 support in McStas / McXtrace

[willend]

TL;DR:

• It is now possible to merge MCPL files from multiple McStas/McXtrace runs to boost statistics (needs weight_mode=1).
• If you use MCPL files from McStas/McXtrace for anything else than feeding them to McStas/McXtrace you may need to rescale weights.
• Going forward, please use MCPL_output with weight_mode=1
• MCPL_input will automatically handle/support both new mode (weight_mode=1) and classical mode (weight_mode=2)
• We will completely transition to weight_mode=1 within 6-12 months.

New stat::sum feature of MCPL

As recently announced by the MCPL team, the latest versions of MCPL includes a new feature for automatic merging of MCPL files with different statistics/normalisation:

This release brings an often-requested feature to MCPL: Statistics in the file header which are automatically combined when files are merged!

The support does not require an update of the on-disk MCPL format, but is handled as a convention for encoding the parameters as key-value pairs in comment fields in the MCPL headers, using a syntax like "stat:sum:<key>:<value>".

The various APIs (C/Python/cmdline) have been updated to provide interfaces for interacting with these data items, and all supported operations (merge/repair/extract/...) have been updated to handle such entries appropriately. For a more complete introduction of the new feature, please refer to the documentation at the https://mctools.github.io/mcpl/format_statsum/ page.

Support for / use of stat::sum in McStas/McXtrace

To fully support this new feature in McStas and McXtrace, the MCPL_output components of both packages have received a new weight_mode input parameter:

• weight_mode=1: record initial ray count in MCPL stat:sum entry and rescale particle weights.
• weight_mode=2: classical (deprecated) mode of outputting particle weights directly.

weight_mode=1 output (the coming default behaviour of MCPL_output)

With weight_mode=1, the events stored in the output MCPL use the stat::sum statistic with a key initial_ray_count, as visible in the MCPL header:

MCPL output summary from tmpdir/1/voutput.mcpl.gz
Opened MCPL file voutput.mcpl.gz:

  Basic info
    Format             : MCPL-3
    No. of particles   : 1000
    Header storage     : 2300 bytes
    Data storage       : 36000 bytes

  Custom meta data
    Source             : "mcstas 3.5.31 Test_MCPL_output"
    Number of comments : 3
          -> comment 0 : "Output by COMPONENT: vout"
          -> comment 1 : "stat:sum:initial_ray_count:                    1000"
          -> comment 2 : "userflags: Neutron Id"
    Number of blobs    : 2
          -> 2013 bytes of data with key "mccode_instr_file"
          -> 48 bytes of data with key "mccode_cmd_line"
...

and particle weights in the table are scaled by the same stat::sum value:

index     pdgcode   ekin[MeV]       x[cm]       y[cm]       z[cm]          ux          uy          uz    time[ms]      weight  userflags
    0        2112  1.1085e-09    -0.44706    -0.25765           0    0.023419   0.0042203     0.99972     0.23609  4.3734e+10 0x00000000
    1        2112  1.5205e-09    -0.39436    -0.33383           0   -0.020454  -0.0035361     0.99978     0.88689  7.8931e+10 0x00000001
    2        2112   1.958e-09    -0.35068    -0.17905           0     0.05127    0.036465     0.99802    0.060739  1.2004e+11 0x00000002
...

weight_mode=2 classic output (the soon-to-be deprecated behaviour of MCPL_output)

With weight_mode=2, there is no stat::sum statistic in the MCPL header:

MCPL output summary from /Users/peterwillendrup/tmp/Test_MCPL_output/Test_MCPL_output_20250617_093316/0/voutput.mcpl.gz
Opened MCPL file voutput.mcpl.gz:

  Basic info
    Format             : MCPL-3
    No. of particles   : 1000
    Header storage     : 2246 bytes
    Data storage       : 36000 bytes

  Custom meta data
    Source             : "mcstas 3.5.31 Test_MCPL_output"
    Number of comments : 2
          -> comment 0 : "Output by COMPONENT: vout"
          -> comment 1 : "userflags: Neutron Id"
    Number of blobs    : 2
          -> 2013 bytes of data with key "mccode_instr_file"
          -> 49 bytes of data with key "mccode_cmd_line"
...

and particle weights in the table are stored directly as they are simulated in McStas/McXtrace:

index     pdgcode   ekin[MeV]       x[cm]       y[cm]       z[cm]          ux          uy          uz    time[ms]      weight  userflags
    0        2112  1.1085e-09    -0.44706    -0.25765           0    0.023419   0.0042203     0.99972     0.23609  4.3734e+07 0x00000000
    1        2112  1.5205e-09    -0.39436    -0.33383           0   -0.020454  -0.0035361     0.99978     0.88689  7.8931e+07 0x00000001
    2        2112   1.958e-09    -0.35068    -0.17905           0     0.05127    0.036465     0.99802    0.060739  1.2004e+08 0x00000002
 ...

Current component default value of -1 explained:

To make it very clear that the classical behaviour will eventually be removed, the default of weight_mode is currently "unspecified" (=-1) which prompts a warning message to take you to this GitHub discussion:

Warning (vout): weight_mode unspecified so defaulting to 2 (classic mode) which will soon become unavailable and tooling might need updating. Please consult https://github.com/mccode-dev/McCode/discussions/2076 for more information.

When will deprecation of weight_mode=2 become effective?

Our plan is to fully transition to weight_mode=1 within 6-12 months. Note that MCPL files written in the legacy format before this will still be readable afterwards, and will not become unusable.

MCPL_input support:

The MCPL_input and MCPL_input_once components both detect presence of the stat::sum field and scale simulated events accordingly, i.e.: no action needed on this end - and future simulations based on old (style) MCPL files will continue to function.

How to update other code reading these files:

If you don't care about the absolute scale of weights in the MCPL files you are reading, you can skip this. Otherwise, you might need to rescale all the weights to get the same meaning of the weights as before. The following Python code shows how this can be done in a way which supports MCPL files from both before and after this change (this needs MCPL release 2.2.0 or later):

def weight_scale( mcplfile ):
    v = mcplfile.stat_sum.get('initial_ray_count',1.0)
    if not v:
        raise RuntimeError('Can not determine weight scale')
    return 1.0 / v

And here is a similar C version:

#include "mcpl.h"
#include <stdio.h>
#include <stdlib.h>

double weight_scale( mcpl_file_t f )
{
  double v = mcpl_hdr_stat_sum( f, "initial_ray_count" );
  if ( v == -2.0 )
    return 1.0;//legacy file
  if ( !(v>0.0) ) {
    fprintf(stderr,"Can not determine weight scale\n");
    exit(1);
  }
  return 1.0 / v;
}

Here is a small example showing it in action on a file with (voutput.mcpl.gz) and without (voutput_legacy.mcpl.gz) the new initial_ray_count feature:

def debug_weights( filename, nmax = 5 ):
    from mcpl import MCPLFile
    f = MCPLFile(filename)
    ws = weight_scale( f )
    print(f'Opened file {filename}'
          f' (determined weight_scale = {ws:g})')
    for i,p in enumerate(f.particles):
        if i==nmax:
            print('   ...')
            break
        print(f'   particle #{i} weights: '
              f'raw={p.weight:g}, rescaled={ws*p.weight:g}')
debug_weights('voutput.mcpl.gz')
debug_weights('voutput_legacy.mcpl.gz')

Opened file voutput.mcpl.gz (determined weight_scale = 0.001)
   particle #0 weights: raw=4.37341e+10, rescaled=4.37341e+07
   particle #1 weights: raw=7.89314e+10, rescaled=7.89314e+07
   particle #2 weights: raw=1.20036e+11, rescaled=1.20036e+08
   particle #3 weights: raw=4.20464e+10, rescaled=4.20464e+07
   particle #4 weights: raw=2.2546e+10, rescaled=2.2546e+07
   ...
Opened file voutput_legacy.mcpl.gz (determined weight_scale = 1)
   particle #0 weights: raw=4.37341e+07, rescaled=4.37341e+07
   particle #1 weights: raw=7.89314e+07, rescaled=7.89314e+07
   particle #2 weights: raw=1.20036e+08, rescaled=1.20036e+08
   particle #3 weights: raw=4.20464e+07, rescaled=4.20464e+07
   particle #4 weights: raw=2.2546e+07, rescaled=2.2546e+07
   ...

Why this page?

Since the implementation touches what the stored weight value in the (MCPL) output files are, we simply want users to be fully aware of the change. :-) In terms of practical use, implications should be very limited for most users.

[willend]

@tkittel, @g5t please comment/correct anything as you see fit. - I am now fixing the warning message in MCPL_output and updating version numbers on the repo toward 3.5.32... (Allowing me to close #2075 and set a new tag after lunch...)