Handling of stat:sum values in MCPL files

.github/resources/conda-setup.yml

@@ -14,7 +14,7 @@ dependencies:
   - bash
   - simple-build-system>=2.0.0
   - ncrystal>=4.1.6
-  - mcpl>=2.0.0
+  - mcpl>=2.2.8
   - matplotlib
   - zlib
   - geant4>=11.4.3,<11.5.0

.github/workflows/condatest.yml

@@ -28,6 +28,11 @@ jobs:
           - { os: ubuntu-latest, CC: gcc-12,   CXX: g++-12,     python: '3.13' }
           - { os: ubuntu-latest, CC: gcc,      CXX: g++,        python: '3.14' }
           - { os: macos-latest,   CC: clang,    CXX: clang++,    python: "3.13" }
+          #Intel (x86_64) macOS, where libc++ compares type information by
+          #address, which makes it sensitive to issues with polymorphic classes
+          #or exceptions shared between Python modules (macos-15-intel is the
+          #last Intel image provided by GitHub, expected until fall 2027):
+          - { os: macos-15-intel, CC: clang,    CXX: clang++,    python: "3.12" }
     name: ${{ matrix.os }}.${{ matrix.CC }}.python-${{ matrix.python }}
     runs-on: ${{ matrix.os }}
     defaults:

doc/source/mcpl.rst

@@ -62,6 +62,18 @@ your :ref:`sim script <sbsimscript>`:
 
 See :ref:`the dedicated documentation <sbparticlegen>` for more options.
 
+If the input file has a ``stat:sum:initial_ray_count`` value, all weights are
+divided by it, so that the weights in the Geant4 simulation are per initial ray
+of the simulation which wrote the file, and do not depend on its number of
+rays. This is for instance the case for files written by the ``MCPL_output``
+component of McStas 3.10 and later (with ``weight_mode=1``, the default), which
+store the weights multiplied by the number of rays. The value of the whole file
+is always used, also when only some of its particles are used (e.g. with
+``skip_events`` or a limited number of events), so the results of jobs which
+each process a part of a file can be added up. Files without the value are used
+with their weights as stored. To always use the weights as stored, set
+``gen.normalise_by_initial_ray_count = False``.
+
 
 MCPL files as simulation output
 -------------------------------
@@ -94,6 +106,51 @@ options, by using the
 module. See an example :sbpkg:`here <G4MCPLTests/scripts/mcplwrite>` for how to
 use it.
 
+The MCPL files written by Geant4 contain ``stat:sum`` values which describe the
+simulation as a whole (independently of how many particles ended up in the
+file), and are added up when files are merged (e.g. with ``mcpltool --merge``,
+or when merging the output of multiple processes):
+
+* ``stat:sum:geant4_nevts``: The number of simulated events.
+* ``stat:sum:initial_geant4_weights``: The sum of the weights of the primary
+  particles of the simulated events, in the same convention as the weights in
+  the file. For instance, the sum of the weights of the particles in the file
+  divided by this value is an efficiency.
+* ``stat:sum:original_initial_ray_count``: The ``stat:sum:initial_ray_count``
+  value of the input file read with MCPLGen (if it has one). This is for
+  reference only: the weights in the file are the weights used in Geant4.
+
+The ``stat:sum:initial_ray_count`` convention of instrument simulations like
+McStas (see above) assumes that each particle in a file is an independent ray.
+This is usually not the case for files written by Geant4, where a single event
+can produce several (correlated) particles. Therefore Geant4 output files do not
+claim the convention by default. If it is appropriate for your use case, the
+``opt_irc`` keyword of ``--mcpl`` (or ``setWriteInitialRayCount(True)`` on an
+``MCPLWriter``) instead records the value of the input file as
+``stat:sum:initial_ray_count``, and writes the weights in the same convention as
+in the input file (i.e. multiplied by the value if MCPLGen divided the weights by
+it). Reading such a file with MCPLGen again then gives the same weights as in the
+original simulation.
+
+When the output of multiple processes is merged, each process writes its share
+of the values describing the whole job (like ``initial_ray_count``), so the
+merged file gets the value.
+
+Custom ``stat:sum`` values describing the whole simulation (for instance the
+duration of a beam pulse) can be added to all MCPL files written by Geant4 from
+the :ref:`sim script <sbsimscript>`, before the simulation is started:
+
+.. code-block:: python
+
+    import G4MCPLPlugins.MCPLWriter as Writer
+    Writer.addStatSum( 'my_beam_time', 3600.0 )
+
+Such values are only added once when the output of multiple processes is
+merged. C++ code can add values computed at the end of the simulation via the
+functions in :sbpkg:`G4MCPL/OutputStatSums.hh<G4MCPL/libinc/OutputStatSums.hh>`,
+which is also how MCPLGen provides the ``stat:sum:initial_ray_count`` of its
+input file.
+
 Modify MCPL files
 -----------------
 
@@ -114,6 +171,12 @@ representative file which can be more easily shared for debugging purposes)::
 
   sb_mcplextra_filterfile infile.mcpl outfile.mcpl -l1000
 
+Any ``stat:sum`` values of the original file (like ``stat:sum:initial_ray_count``
+from McStas) are kept when filtering particles, since the new file is simply a
+narrower view of the same simulation. But when only a part of the particles is
+used due to the limit (``-l``), the values are marked as unavailable, since the
+file was then split rather than filtered.
+
 Run ``sb_mcplextra_filterfile --help`` for full usage instructions. Note that
 there also exists a ``sb_mcplextra_filterview`` command, which merely prints the
 selected particles, rather than producing a new file.

pyproject.toml

@@ -28,13 +28,13 @@ dependencies = [
 recommended = [
     'numpy',
     'ncrystal>=4.2.10',
-    'mcpl>=2.2.4',
+    'mcpl>=2.2.8',
     'matplotlib'
 ]
 all = [
     'numpy',
     'ncrystal[all]>=4.2.10',
-    'mcpl>=2.2.4',
+    'mcpl>=2.2.8',
     'matplotlib',
     'scipy',
     'pandas',

src/simplebuild_dgcode/data/pkgs/G4/G4Launcher/python/_launcher.py

@@ -232,6 +232,10 @@ def _display_mcpl_help():
     print("   opt_dp        : Store double-precision numbers in output.")
     print("   opt_pol       : Store polarisation info in output.")
     print("   opt_userflags : Store special user flags in output.")
+    print("   opt_irc       : Store stat:sum:initial_ray_count of the MCPL input")
+    print("                   file, with weights in the same convention as in")
+    print("                   that file (only if each particle can be treated")
+    print("                   as an independent ray).")
     print("   grabsrc       : Short-cut to grab just particles from the particle")
     print('                   generator. Implies kill_never and all_vols plus')
     print('                   prepends "trk.is_primary && step.stepnbr==1" to')
@@ -352,7 +356,7 @@ def _configure_mcplwriter(launcher,cfgstr,parseerror):
     ks_default='NEVER' if 'grabsrc' in cfgs else 'ALWAYS'#default to avoid unintentional double-counting
     on_exit_default = False
     ks, on_exit = [], []
-    opt_dp,opt_pol,opt_userflags = False,False,False
+    opt_dp,opt_pol,opt_userflags,opt_irc = False,False,False,False
     for cfg in cfgs:
         cfg = cfg.strip()
         if not cfg:
@@ -370,6 +374,7 @@ def _configure_mcplwriter(launcher,cfgstr,parseerror):
         elif cfg=='opt_dp': opt_dp = True
         elif cfg=='opt_pol': opt_pol = True
         elif cfg=='opt_userflags': opt_userflags = True
+        elif cfg=='opt_irc': opt_irc = True
         elif cfg=='all_vols': mw.addVolume('*')
         elif cfg == 'where' or cfg == 'to' or cfg.startswith(('kill_', 'on_vol_', 'opt_')) or cfg == 'withflags':
             parseerror('Suspicious volume name "%s" indicates malformed option for --mcpl'%cfg)
@@ -380,6 +385,7 @@ def _configure_mcplwriter(launcher,cfgstr,parseerror):
     mw.setWriteDoublePrecision(opt_dp)
     mw.setWritePolarisation(opt_pol)
     mw.setWriteUserFlags(opt_userflags)
+    mw.setWriteInitialRayCount(opt_irc)
     mw.setWriteOnVolExit(on_exit[0] if on_exit else on_exit_default)
     mw.setKillStrategy(ks[0] if ks else ks_default)
     if mcpl_fltr:

src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/OutputStatSums.hh

@@ -0,0 +1,58 @@
+#ifndef G4MCPL_OutputStatSums_hh
+#define G4MCPL_OutputStatSums_hh
+
+#include <string>
+#include <vector>
+#include <functional>
+
+//Global (per process) registry of information for the MCPL files written by
+//G4MCPLPlugins/MCPLWriter:
+//
+//  * stat:sum values to add to all its output files (e.g. from user code).
+//  * the stat:sum:initial_ray_count of an MCPL input file (set by
+//    G4MCPLPlugins/MCPLGen), which MCPLWriter records in its output files.
+//
+//Since MCPL requires the same header entries in all files to merge them (e.g.
+//the files of the processes of multi-process jobs), the information must be
+//provided before the output files are created, i.e. before the simulation
+//starts (for instance when the parameters of a generator are validated, or in
+//a pre-init hook).
+
+namespace G4MCPL {
+
+  //How a value adds up when the output files of several processes are merged
+  //(MCPL adds up stat:sum values when merging files):
+  enum class StatSumMode {
+    PerProcess,//Each process writes its own value (e.g. a count of events).
+    WholeJob   //The value describes the whole job, so each of the N processes
+               //writes the value divided by N (the merged file then gets the
+               //value). An unavailable value (-1) is written as -1.
+  };
+
+  //Add a stat:sum entry. The function is evaluated when each output file is
+  //closed, and must return a non-negative value, or -1 if the value is not
+  //available. The key must be a valid MCPL stat:sum key, which is not already
+  //added (or used by MCPLWriter itself):
+  void addOutputStatSum( const std::string& key,
+                         std::function<double()> value,
+                         StatSumMode mode );
+
+  struct OutputStatSum {
+    std::string key;
+    std::function<double()> value;
+    StatSumMode mode;
+  };
+  const std::vector<OutputStatSum>& outputStatSums();
+
+  //The stat:sum:initial_ray_count value of an MCPL input file (positive, or -1
+  //if not available), and the factor which was applied to the weights of the
+  //input file when generating the particles (1/initial_ray_count when
+  //normalising, otherwise 1). Set by MCPLGen:
+  void setInputInitialRayCount( double initial_ray_count, double weight_factor );
+  bool hasInputInitialRayCount();
+  double inputInitialRayCount();
+  double inputWeightFactor();
+
+}
+
+#endif

src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc

@@ -0,0 +1,66 @@
+#include "G4MCPL/OutputStatSums.hh"
+#include <stdexcept>
+#include <cmath>
+
+namespace G4MCPL {
+  namespace {
+    std::vector<OutputStatSum>& registry()
+    {
+      static std::vector<OutputStatSum> s_registry;
+      return s_registry;
+    }
+    struct InputInfo {
+      bool present = false;
+      double initial_ray_count = -1.0;
+      double weight_factor = 1.0;
+    };
+    InputInfo& inputInfo()
+    {
+      static InputInfo s_info;
+      return s_info;
+    }
+  }
+}
+
+void G4MCPL::addOutputStatSum( const std::string& key,
+                               std::function<double()> value,
+                               StatSumMode mode )
+{
+  //(the remaining requirements on keys are checked by MCPL when the files are
+  //created):
+  if ( key.empty() || key.size() > 64 || key.find(':') != std::string::npos )
+    throw std::runtime_error("G4MCPL: Invalid stat:sum key (must be non-empty, at"
+                             " most 64 characters, and not contain colons): \""+key+"\"");
+  if ( key == "geant4_nevts" || key == "initial_geant4_weights"
+       || key == "initial_ray_count" || key == "original_initial_ray_count" )
+    throw std::runtime_error("G4MCPL: The stat:sum key \""+key+"\" is reserved for MCPLWriter");
+  for ( auto& e : registry() )
+    if ( e.key == key )
+      throw std::runtime_error("G4MCPL: The stat:sum key \""+key+"\" was already added");
+  if ( !value )
+    throw std::runtime_error("G4MCPL: No function given for stat:sum key \""+key+"\"");
+  registry().push_back( OutputStatSum{ key, std::move(value), mode } );
+}
+
+const std::vector<G4MCPL::OutputStatSum>& G4MCPL::outputStatSums()
+{
+  return registry();
+}
+
+void G4MCPL::setInputInitialRayCount( double initial_ray_count, double weight_factor )
+{
+  if ( !( weight_factor > 0.0 ) || !std::isfinite(weight_factor) )
+    throw std::runtime_error("G4MCPL: Weight factor must be positive and finite");
+  if ( !( initial_ray_count > 0.0 ) || !std::isfinite(initial_ray_count) )
+    initial_ray_count = -1.0;//not available
+  auto& info = inputInfo();
+  if ( info.present )
+    throw std::runtime_error("G4MCPL: The initial_ray_count of an input file was already set");
+  info.present = true;
+  info.initial_ray_count = initial_ray_count;
+  info.weight_factor = weight_factor;
+}
+
+bool G4MCPL::hasInputInitialRayCount() { return inputInfo().present; }
+double G4MCPL::inputInitialRayCount() { return inputInfo().initial_ray_count; }
+double G4MCPL::inputWeightFactor() { return inputInfo().weight_factor; }