diff --git a/.github/resources/conda-setup.yml b/.github/resources/conda-setup.yml index c2d6a38d..16cbac70 100644 --- a/.github/resources/conda-setup.yml +++ b/.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 diff --git a/.github/workflows/condatest.yml b/.github/workflows/condatest.yml index 373b6d25..714012a3 100644 --- a/.github/workflows/condatest.yml +++ b/.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: diff --git a/doc/source/mcpl.rst b/doc/source/mcpl.rst index 1ebe9273..ac00af4c 100644 --- a/doc/source/mcpl.rst +++ b/doc/source/mcpl.rst @@ -62,6 +62,18 @@ your :ref:`sim script `: See :ref:`the dedicated documentation ` 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 ` 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 `, 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`, +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. diff --git a/pyproject.toml b/pyproject.toml index f1db4b63..a36c8ad8 100644 --- a/pyproject.toml +++ b/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', diff --git a/src/simplebuild_dgcode/data/pkgs/G4/G4Launcher/python/_launcher.py b/src/simplebuild_dgcode/data/pkgs/G4/G4Launcher/python/_launcher.py index 95a7110d..6f57a55b 100644 --- a/src/simplebuild_dgcode/data/pkgs/G4/G4Launcher/python/_launcher.py +++ b/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: diff --git a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/OutputStatSums.hh b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/OutputStatSums.hh new file mode 100644 index 00000000..fb2f6375 --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/OutputStatSums.hh @@ -0,0 +1,58 @@ +#ifndef G4MCPL_OutputStatSums_hh +#define G4MCPL_OutputStatSums_hh + +#include +#include +#include + +//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 value, + StatSumMode mode ); + + struct OutputStatSum { + std::string key; + std::function value; + StatSumMode mode; + }; + const std::vector& 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 diff --git a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc new file mode 100644 index 00000000..e5368a04 --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc @@ -0,0 +1,66 @@ +#include "G4MCPL/OutputStatSums.hh" +#include +#include + +namespace G4MCPL { + namespace { + std::vector& registry() + { + static std::vector 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 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::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; } diff --git a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc index 3bd7ea99..671946b7 100644 --- a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc @@ -1,9 +1,13 @@ #include "G4MCPL/G4MCPLUserFlags.hh" +#include "G4MCPL/OutputStatSums.hh" #include "G4Interfaces/ParticleGenToPy.hh" #include "G4Interfaces/ParticleGenBase.hh" #include "G4Interfaces/FrameworkGlobals.hh" #include "MCPLExprParser/MCPLASTBuilder.hh" #include "mcpl.h" +#if !defined(MCPL_VERSION) || MCPL_VERSION < 2002008 +# error "MCPL version 2.2.8 or later is required" +#endif #include "Core/FindData.hh" #include "Units/Units.hh" #include "Core/String.hh" @@ -12,6 +16,7 @@ #include "CLHEP/Geometry/Vector3D.h" #include "CLHEP/Geometry/Transform3D.h" #include +#include #include #include #include "G4ParticleTable.hh" @@ -33,6 +38,16 @@ class MCPLGen : public G4Interfaces::ParticleGenBase { void skip_forward(unsigned nskip); void delayed_init(); void setPDG(int); + void init_weight_scale(); + void register_output_stat_sums(); + std::string resolved_input_file() const; + //The stat:sum:initial_ray_count value of the input file (-2 if not present, + //-1 if not available), the factor applied to its weights (1/initial_ray_count + //when normalising, otherwise 1), and whether all of its particles were used + //(in this process, when multiprocessing): + double inputInitialRayCount() const; + double inputWeightFactor() const; + mcpl_file_t m_mcplfile; const mcpl_particle_t * m_p;//particle to be used for next generation bool m_unfiltered; @@ -48,6 +63,9 @@ class MCPLGen : public G4Interfaces::ParticleGenBase { HepGeom::Vector3D m_translation; HepGeom::Transform3D m_rotation; bool m_allow_zero_weight; + double m_weight_scale;//factor applied to all weights in the file + mutable double m_input_irc;//cache for inputInitialRayCount() (NaN if not loaded) + bool m_registered_output_stat_sums; }; MCPLGen::MCPLGen() @@ -58,7 +76,10 @@ MCPLGen::MCPLGen() m_nprocs(0), m_lastpdgcode(0), m_has_polarisation(false), - m_allow_zero_weight(false) + m_allow_zero_weight(false), + m_weight_scale(1.0), + m_input_irc(std::nan("")), + m_registered_output_stat_sums(false) { m_mcplfile.internal = 0; @@ -83,6 +104,14 @@ MCPLGen::MCPLGen() //If not set, input events with weight=0 triggers an error. addParameterBoolean("allow_zero_weight",false); + //Divide the weights by stat:sum:initial_ray_count if the file has it + //(e.g. files from McStas MCPL_output with weight_mode=1, the default since + //McStas 3.10), so that the weights are per initial ray of the simulation + //that wrote the file and do not depend on its ray count. The value of the + //whole file is used, also when only some of its particles are used (e.g. + //with skip_events, or a limited number of events): + addParameterBoolean("normalise_by_initial_ray_count",true); + } MCPLGen::~MCPLGen() @@ -119,9 +148,34 @@ bool MCPLGen::validateParameters() m_unfiltered = ( m_eval_filter.isConstant() && m_eval_filter() ); + //The input file must be found here (rather than when it is opened for reading + //the particles), since its header is needed for registering the values of + //the output files: + if (Core::findData(fn).empty()) { + printf("MCPLGen ERROR - File specified in input_file parameter not found : \"%s\"\n",fn.c_str()); + return false; + } + register_output_stat_sums(); + return true; } +void MCPLGen::register_output_stat_sums() +{ + //Provide the stat:sum:initial_ray_count of the input file (and the factor + //applied to its weights) to MCPLWriter, which records it in its output files + //(see G4MCPL/OutputStatSums.hh). This happens when the parameters are + //validated, which is before the simulation starts (and before any fork when + //multiprocessing): + if (m_registered_output_stat_sums) + return; + m_registered_output_stat_sums = true; + const double nrays = inputInitialRayCount(); + if (nrays == -2.0) + return;//not present + G4MCPL::setInputInitialRayCount( nrays, inputWeightFactor() ); +} + void MCPLGen::init() { //In case of multiprocessing we haven't forked yet and each process needs it's @@ -157,20 +211,51 @@ void MCPLGen::skip_forward(unsigned nskip) { } } -void MCPLGen::delayed_init() +std::string MCPLGen::resolved_input_file() const { - //opening files should happen after fork() when multi-processing. std::string fn = getParameterString("input_file"); std::string fn_resolved = Core::findData(fn); if (fn_resolved.empty()) { printf("MCPLGen ERROR - File specified in input_file parameter not found : \"%s\"\n",fn.c_str()); throw std::runtime_error("MCPLGen: File specified in input_file parameter not found"); } + return fn_resolved; +} + +double MCPLGen::inputInitialRayCount() const +{ + if (std::isnan(m_input_irc)) { + if (m_mcplfile.internal) { + m_input_irc = mcpl_hdr_stat_sum(m_mcplfile,"initial_ray_count"); + } else { + //Only read the header (the file is opened for reading the particles + //later, after any fork() when multiprocessing): + mcpl_file_t f = mcpl_open_file(resolved_input_file().c_str()); + m_input_irc = mcpl_hdr_stat_sum(f,"initial_ray_count"); + mcpl_close_file(f); + } + } + return m_input_irc; +} - m_mcplfile = mcpl_open_file(fn_resolved.c_str()); +double MCPLGen::inputWeightFactor() const +{ + const double nrays = inputInitialRayCount(); + if ( getParameterBoolean("normalise_by_initial_ray_count") && nrays > 0.0 && std::isfinite(nrays) ) + return 1.0 / nrays; + return 1.0; +} + + +void MCPLGen::delayed_init() +{ + //opening files should happen after fork() when multi-processing. + m_mcplfile = mcpl_open_file(resolved_input_file().c_str()); m_has_polarisation = mcpl_hdr_has_polarisation(m_mcplfile); + init_weight_scale(); + unsigned nskip = getParameterInt("skip_events"); if (nskip) skip_forward(nskip); @@ -191,6 +276,43 @@ void MCPLGen::delayed_init() } } +void MCPLGen::init_weight_scale() +{ + m_weight_scale = 1.0; + const double nrays = inputInitialRayCount(); + //Only print in one process when multiprocessing: + const bool doprint = FrameworkGlobals::isParent(); + if (!getParameterBoolean("normalise_by_initial_ray_count")) { + if (doprint && nrays != -2.0) { + printf("MCPLGen: Using the weights in the file as stored (normalise_by_initial_ray_count" + " is disabled, stat:sum:initial_ray_count is %g)\n",nrays); + fflush(stdout); + } + return; + } + if (nrays == -2.0) + return;//not present (e.g. files from older McStas versions): use weights as stored + if ( !(nrays > 0.0) || !std::isfinite(nrays) ) { + if (nrays == -1.0) + printf("MCPLGen ERROR - The input file has a stat:sum:initial_ray_count value which is not" + " available (for instance due to an interrupted run of the application writing the file," + " or because the file contains only a part of the particles of another file), so the weights" + " can not be normalised. Set normalise_by_initial_ray_count=false to use the weights as" + " stored in the file.\n"); + else + printf("MCPLGen ERROR - The input file has an invalid stat:sum:initial_ray_count value (%g)," + " so the weights can not be normalised. Set normalise_by_initial_ray_count=false to use" + " the weights as stored in the file.\n",nrays); + throw std::runtime_error("MCPLGen: Can not normalise weights by stat:sum:initial_ray_count"); + } + m_weight_scale = inputWeightFactor(); + if (doprint) { + printf("MCPLGen: Dividing all weights by stat:sum:initial_ray_count = %g" + " (disable with normalise_by_initial_ray_count=false)\n",nrays); + fflush(stdout); + } +} + void MCPLGen::setPDG(int p) { assert(m_lastpdgcode!=p); @@ -273,11 +395,13 @@ void MCPLGen::gen(G4Event* evt) assert(evt->GetNumberOfPrimaryVertex()==0); const int ivertex = 0; m_gun->GeneratePrimaryVertex(evt); - if (m_p->weight!=1.0) { - if ( !m_p->weight && !m_allow_zero_weight ) { + //(the weights are normalised before checking if they are 1): + const double weight = m_p->weight * m_weight_scale; + if (weight!=1.0) { + if ( !weight && !m_allow_zero_weight ) { throw std::runtime_error("MCPLGen: File contains particles with weight=0 and allow_zero_weight was not enabled."); } - evt->GetPrimaryVertex(ivertex)->SetWeight(m_p->weight); + evt->GetPrimaryVertex(ivertex)->SetWeight(weight); } //User flags require a bit special treatment: diff --git a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLWriter/mod.cc b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLWriter/mod.cc index 2b813b02..af7dc4f1 100644 --- a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLWriter/mod.cc +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLWriter/mod.cc @@ -2,18 +2,27 @@ #include "G4Interfaces/ComponentHandles.hh" #include "G4MCPL/G4MCPLUserFlags.hh" +#include "G4MCPL/OutputStatSums.hh" #include "G4Interfaces/FrameworkGlobals.hh" #include "G4Interfaces/GeoConstructBase.hh" #include "G4Interfaces/ParticleGenBase.hh" #include "G4ExprParser/G4SteppingASTBuilder.hh" #include "mcpl.h" +#if !defined(MCPL_VERSION) || MCPL_VERSION < 2002008 +# error "MCPL version 2.2.8 or later is required" +#endif #include "G4VSensitiveDetector.hh" +#include "G4Event.hh" +#include "G4PrimaryVertex.hh" +#include "G4PrimaryParticle.hh" #include "G4Step.hh" #include "G4SDManager.hh" #include "G4RegionStore.hh" #include +#include +#include #include #include #include @@ -34,6 +43,49 @@ namespace G4MCPLWriter { MCPLSensitiveDetector * m_sd; }; + //Statistics of the simulated events, which are written to the output files + //as stat:sum values (they are added up when files are merged): + // + // geant4_nevts : The number of simulated events. + // initial_geant4_weights : The sum of the weights of the primary particles + // of the simulated events. + // + //They are collected via a post-generation callback on the particle + //generator (so they are independent of which particles reach the volumes). + class EventStats : public G4Interfaces::PostGenCallBack { + public: + EventStats( const G4Interfaces::ParticleGenBase* gen ) : m_gen(gen) {} + void postGen( G4Event* evt ) override + { + const int nvertices = evt->GetNumberOfPrimaryVertex(); + if ( !nvertices && m_gen->reachedLimit() ) + return;//the generator ran out of events, so this event is aborted + ++m_nevts; + for ( int i = 0; i < nvertices; ++i ) { + const G4PrimaryVertex* vtx = evt->GetPrimaryVertex(i); + for ( const G4PrimaryParticle* pp = vtx->GetPrimary(); pp; pp = pp->GetNext() ) + m_sumweights += sumWeights( pp, vtx->GetWeight() ); + } + } + std::uint64_t nevts() const { return m_nevts; } + double sumWeights() const { return m_sumweights; } + private: + static double sumWeights( const G4PrimaryParticle* pp, double vtxweight ) + { + //The weight of the track of a primary particle is the product of the + //weights of the vertex and the particle (also for pre-assigned + //daughters): + double w = vtxweight * pp->GetWeight(); + for ( const G4PrimaryParticle* d = pp->GetDaughter(); d; d = d->GetNext() ) + w += sumWeights( d, vtxweight ); + return w; + } + const G4Interfaces::ParticleGenBase* m_gen; + std::uint64_t m_nevts = 0; + double m_sumweights = 0.0; + }; + typedef std::shared_ptr EventStatsPtr; + enum KillStrategy { KS_NEVER, KS_ALWAYS, KS_FILTERED }; KillStrategy ks_from_string(const std::string& s) { @@ -68,12 +120,14 @@ namespace G4MCPLWriter { bool writepolarisation, bool writeuserflags, double universalweight, + bool writeinitialraycount, bool writeonvolexit, KillStrategy killstrategy, BuilderPtr builder, EP::Evaluator filter_evaluator, EP::Evaluator flag_evaluator, - std::vector> comments_and_blobs) + std::vector> comments_and_blobs, + EventStatsPtr evtstats) : G4VSensitiveDetector(name), m_expr_builder(builder), m_eval_filter(filter_evaluator), @@ -89,7 +143,10 @@ namespace G4MCPLWriter { m_initialised(false), m_closed(false), m_om(0), - m_comments_and_blobs(comments_and_blobs) + m_comments_and_blobs(comments_and_blobs), + m_evtstats(evtstats), + m_opt_writeinitialraycount(writeinitialraycount), + m_weight_factor(1.0) { std::memset(&m_p,0,sizeof(m_p)); } @@ -125,6 +182,21 @@ namespace G4MCPLWriter { return; m_closed = true; assert(m_f.internal); + //Final values of the stat:sum entries reserved in initmcpl(): + if (m_evtstats) { + mcpl_hdr_add_stat_sum(m_f,"geant4_nevts",(double)m_evtstats->nevts()); + //(in the same convention as the weights in the file): + mcpl_hdr_add_stat_sum(m_f,"initial_geant4_weights",m_evtstats->sumWeights()*m_weight_factor); + } + if (!m_irc_key.empty()) + mcpl_hdr_add_stat_sum(m_f,m_irc_key.c_str(),wholeJobValue(G4MCPL::inputInitialRayCount())); + for ( auto& e : m_statsums ) { + const double value = e.value(); + mcpl_hdr_add_stat_sum(m_f,e.key.c_str(), + ( e.mode == G4MCPL::StatSumMode::WholeJob + ? wholeJobValue(value) + : ( value >= 0.0 && std::isfinite(value) ? value : -1.0 ) )); + } if (!FrameworkGlobals::isForked()) mcpl_closeandgzip_outfile(m_f); else @@ -195,7 +267,7 @@ namespace G4MCPLWriter { } m_p.time = pt->GetGlobalTime() / CLHEP::millisecond; - m_p.weight = pt->GetWeight(); + m_p.weight = pt->GetWeight() * m_weight_factor; assert(CLHEP::MeV==1.0); m_p.ekin = pt->GetKineticEnergy();//already in MeV const G4ThreeVector& pos = pt->GetPosition(); @@ -268,10 +340,40 @@ namespace G4MCPLWriter { mcpl_hdr_add_data(m_f,it->first.c_str(),it->second.size(),it->second.c_str()); } + if (m_evtstats) { + //Reserve space for the stat:sum values, which are set when the file is + //closed (see ensure_close_file()): + mcpl_hdr_add_stat_sum(m_f,"geant4_nevts",-1.0); + mcpl_hdr_add_stat_sum(m_f,"initial_geant4_weights",-1.0); + } + + //The stat:sum:initial_ray_count of an input file read by MCPLGen. By + //default it is recorded as original_initial_ray_count, with the weights + //as used in Geant4. Only if requested, it is recorded as + //initial_ray_count, with the weights in the same convention as in the + //input file (i.e. multiplied by the value if MCPLGen divided the weights + //by it), so that reading the file again gives the same weights. That + //convention assumes that each particle is an independent ray, which is + //for the user to decide. The space for the values is reserved here, and + //the values are set when the file is closed: + m_irc_key.clear(); + m_weight_factor = 1.0; + if ( G4MCPL::hasInputInitialRayCount() ) { + m_irc_key = ( m_opt_writeinitialraycount ? "initial_ray_count" : "original_initial_ray_count" ); + if ( m_opt_writeinitialraycount ) + m_weight_factor = 1.0 / G4MCPL::inputWeightFactor(); + mcpl_hdr_add_stat_sum(m_f,m_irc_key.c_str(),-1.0); + } + + //Values added by other components or user code: + m_statsums = G4MCPL::outputStatSums(); + for ( auto& e : m_statsums ) + mcpl_hdr_add_stat_sum(m_f,e.key.c_str(),-1.0); + if (m_opt_writedoubleprec) mcpl_enable_doubleprec(m_f); if (m_opt_writepolarisation) mcpl_enable_polarisation(m_f); if (m_opt_writeuserflags) mcpl_enable_userflags(m_f); - if (m_opt_universalweight) mcpl_enable_universal_weight(m_f,m_opt_universalweight); + if (m_opt_universalweight) mcpl_enable_universal_weight(m_f,m_opt_universalweight*m_weight_factor); if (FrameworkGlobals::isForked()&&FrameworkGlobals::isParent()) m_om = new MCPLOutputMerger(this);//used by mergeIfNeeded() @@ -290,6 +392,25 @@ namespace G4MCPLWriter { mcpl_particle_t m_p; MCPLOutputMerger * m_om; std::vector> m_comments_and_blobs; + EventStatsPtr m_evtstats; + std::vector m_statsums; + bool m_opt_writeinitialraycount; + std::string m_irc_key; + double m_weight_factor; + + static double wholeJobValue( double value ) + { + //A value describing the whole job. When multiprocessing, each of the N + //processes writes the value divided by N, so the merged file gets the + //value (MCPL adds up stat:sum values when merging files). + // + //NB: The value would preferably be set only in the merged file, but MCPL + //does not currently support updating stat:sum values of existing files + //(see mctools/mcpl#254 and mctools/simplebuild-dgcode#63). + if ( !( value >= 0.0 ) || !std::isfinite(value) ) + return -1.0;//not available + return value / std::max( 1, FrameworkGlobals::nProcs() ); + } }; void MCPLOutputMerger::merge() { @@ -321,6 +442,7 @@ namespace G4MCPLWriter { bool m_opt_writeonvolexit; KillStrategy m_opt_killstrategy; double m_opt_universalweight; + bool m_opt_writeinitialraycount; BuilderPtr m_expr_builder; EP::Evaluator m_eval_filter; EP::Evaluator m_eval_flags; @@ -332,6 +454,7 @@ namespace G4MCPLWriter { void setWriteUserFlags(bool b) { m_opt_writeuserflags = b; } void setWriteOnVolExit(bool b) { m_opt_writeonvolexit = b; } void setUniversalWeight(double w) { m_opt_universalweight = w; } + void setWriteInitialRayCount(bool b) { m_opt_writeinitialraycount = b; } void setFilter(const char * expr) { m_expr_filter = expr; } void setUserFlags(const char * expr) { m_expr_flags = expr; } void setKillStrategy(KillStrategy ks) { m_opt_killstrategy = ks; } @@ -345,6 +468,7 @@ namespace G4MCPLWriter { m_opt_writeonvolexit(false), m_opt_killstrategy(KS_NEVER), m_opt_universalweight(0.0), + m_opt_writeinitialraycount(false), m_expr_filter("true"), m_vols_all(false), m_sd(0) @@ -370,10 +494,20 @@ namespace G4MCPLWriter { m_geo = geo_handle ? &G4Interfaces::getComponentAs( geo_handle, "geometry" ) : nullptr; m_gen = gen_handle ? &G4Interfaces::getComponentAs( gen_handle, "particle generator" ) : nullptr; + if ( m_gen && !m_evtstats ) { + m_evtstats = std::make_shared( m_gen ); + m_gen->installPostGenCallBack( m_evtstats ); + } if (m_vols.empty()&&!m_vols_all) throw std::runtime_error("MCPLWriter: No volumes added!"); + if ( m_opt_writeinitialraycount && !( G4MCPL::hasInputInitialRayCount() + && G4MCPL::inputInitialRayCount() > 0.0 ) ) + throw std::runtime_error("MCPLWriter: Writing stat:sum:initial_ray_count was requested, but there" + " is no MCPL input file (read with MCPLGen) with a valid" + " stat:sum:initial_ray_count value"); + //Setup filter: if (m_eval_filter.arg()|| m_expr_builder) @@ -448,7 +582,8 @@ namespace G4MCPLWriter { bool m_vols_all; MCPLSensitiveDetector * m_sd; const G4Interfaces::GeoConstructBase* m_geo = nullptr; - const G4Interfaces::ParticleGenBase* m_gen = nullptr; + G4Interfaces::ParticleGenBase* m_gen = nullptr; + EventStatsPtr m_evtstats; std::vector> createCommentsAndBlobs() const { std::vector> out; @@ -536,6 +671,11 @@ namespace G4MCPLWriter { tmp << "ESS/dgcode generator module : " << gen_name; out.emplace_back("comment",tmp.str()); } + if (gen) + out.emplace_back("comment","MCPLWriter stat:sum values : geant4_nevts is the number of simulated" + " events, initial_geant4_weights is the sum of the weights of their primary particles" + " (in the convention of the weights in this file), original_initial_ray_count is the" + " initial_ray_count of the input file (if any)"); if (geo) { char * geodata(0); @@ -583,12 +723,14 @@ namespace G4MCPLWriter { m_opt_writepolarisation, m_opt_writeuserflags, m_opt_universalweight, + m_opt_writeinitialraycount, m_opt_writeonvolexit, m_opt_killstrategy, m_expr_builder, m_eval_filter, m_eval_flags, - createCommentsAndBlobs()); + createCommentsAndBlobs(), + m_evtstats); G4SDManager::GetSDMpointer()->AddNewDetector(m_sd); } auto currentsd = lv->GetSensitiveDetector(); @@ -647,7 +789,13 @@ SBLDPY_MODULE( _MCPLWriter ) .def("setWritePolarisation",[]( PyMCPLWriter& o, bool b ) { o.w().setWritePolarisation( b ); }, { "b" } ) .def("setWriteUserFlags",[]( PyMCPLWriter& o, bool b ) { o.w().setWriteUserFlags( b ); }, { "b" } ) .def("setUniversalWeight",[]( PyMCPLWriter& o, double w ) { o.w().setUniversalWeight( w ); }, { "w" } ) + .def("setWriteInitialRayCount",[]( PyMCPLWriter& o, bool b ) { o.w().setWriteInitialRayCount( b ); }, { "b" } ) .def("setKillStrategy",[]( PyMCPLWriter& o, const std::string& ks ) { o.w().setKillStrategyFromString( ks.c_str() ); }, { "ks" } ) .def("inithook",[]( PyMCPLWriter& o, std::int64_t geo, std::int64_t gen ) { o.w().inithook( geo, gen ); }, { "geo_handle", "gen_handle" } ) .def("postmp",[]( PyMCPLWriter& o ) { o.w().postmp(); } ); + //Add a constant stat:sum value describing the whole job to all output files + //(see G4MCPL/OutputStatSums.hh): + m.def("addStatSum",[]( const std::string& key, double value ) { + G4MCPL::addOutputStatSum( key, [value](){ return value; }, G4MCPL::StatSumMode::WholeJob ); + }, { "key", "value" } ); } diff --git a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/python/MCPLWriter.py b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/python/MCPLWriter.py index ead7478a..c950e53d 100644 --- a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/python/MCPLWriter.py +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/python/MCPLWriter.py @@ -1,4 +1,13 @@ -__all__ = [ 'MCPLWriter' ] +__all__ = [ 'MCPLWriter', 'addStatSum' ] + +def addStatSum( key, value ): + """Add a stat:sum entry with a constant value describing the whole + simulation (e.g. the duration of a beam pulse), to all MCPL files written by + MCPLWriter (including those written due to the --mcpl option). Must be + called before the simulation starts. Use -1 to mark the value as not + available.""" + from ._MCPLWriter import addStatSum as _addStatSum + _addStatSum( key, float(value) ) class MCPLWriter: """Capture particles in a Geant4 simulation and write them to an MCPL file. @@ -10,7 +19,7 @@ def __init__( self, filename ): self._w = _MCPLWriter( filename ) for n in ( 'addVolume', 'setFilter', 'setUserFlags', 'setWriteOnVolExit', 'setWriteDoublePrecision', 'setWritePolarisation', 'setWriteUserFlags', - 'setUniversalWeight', 'setKillStrategy' ): + 'setUniversalWeight', 'setKillStrategy', 'setWriteInitialRayCount' ): setattr( self, n, getattr( self._w, n ) ) def inithook( self ): diff --git a/src/simplebuild_dgcode/data/pkgs/MCPL/MCPLExtra/app_filterfile/main.cc b/src/simplebuild_dgcode/data/pkgs/MCPL/MCPLExtra/app_filterfile/main.cc index 12405d90..4d94fd9a 100644 --- a/src/simplebuild_dgcode/data/pkgs/MCPL/MCPLExtra/app_filterfile/main.cc +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/MCPLExtra/app_filterfile/main.cc @@ -168,6 +168,7 @@ int main(int argc, char** argv) { double norig = mcpl_hdr_nparticles(fi); int progress = -1; const mcpl_particle_t* p; + bool truncated = false; printf("Start processing particle data.\n"); while ( ( p=mcpl_read(fi) ) ) { ++posp1; @@ -183,11 +184,26 @@ int main(int argc, char** argv) { continue; mcpl_add_particle(fo,p); ++used; - if ( nparticles_limit && used==nparticles_limit ) + if ( nparticles_limit && used==nparticles_limit ) { + truncated = ( mcpl_currentposition(fi) < mcpl_hdr_nparticles(fi) ); break; + } } printf("Done processing particle data.\n"); + //Any stat:sum values (like stat:sum:initial_ray_count from McStas) remain + //valid when filtering the particles, but not when only a part of the file is + //used due to the limit (the file was split, rather than filtered): + bool has_stat_sums = false; + for ( unsigned i = 0; i < mcpl_hdr_ncomments(fi); ++i ) + if ( std::strncmp( mcpl_hdr_comment(fi,i), "stat:sum:", 9 ) == 0 ) + has_stat_sums = true; + if ( truncated && has_stat_sums ) { + printf("Marking stat:sum values as unavailable, since the limit on the number" + " of particles means that only a part of the input file was used.\n"); + mcpl_hdr_scale_stat_sums(fo,-1.0); + } + //close: if (do_gzip) mcpl_closeandgzip_outfile(fo); diff --git a/src/simplebuild_dgcode/data/pkgs/Utils/ExprParser/libinc/Exception.hh b/src/simplebuild_dgcode/data/pkgs/Utils/ExprParser/libinc/Exception.hh index db5c6829..d4033ad8 100644 --- a/src/simplebuild_dgcode/data/pkgs/Utils/ExprParser/libinc/Exception.hh +++ b/src/simplebuild_dgcode/data/pkgs/Utils/ExprParser/libinc/Exception.hh @@ -6,12 +6,20 @@ namespace ExprParser { - //Base class for all ExprParser exceptions: + //Base class for all ExprParser exceptions. + // + //NB: The destructors of all exception classes are defined in Exception.cc + //(rather than inline), so that their type information is only present in the + //ExprParser library. Otherwise each module using them (e.g. Python modules, + //which are loaded with RTLD_LOCAL) would get its own copy, and on some + //platforms (e.g. Intel macOS, where libc++ compares type information by + //address) an exception thrown in one module could not be caught by type in + //another module. class Exception { public: Exception() {} - virtual ~Exception() throw() {} + virtual ~Exception() throw(); virtual const char * epType() const = 0; virtual const char * epWhat() const = 0; }; @@ -24,6 +32,7 @@ namespace ExprParser { class InputError : public Exception { public: using Exception::Exception; + virtual ~InputError() throw(); }; //Exception thrown if input expressions contains errors (for instance syntax @@ -32,6 +41,7 @@ namespace ExprParser { class ParseError : public std::runtime_error, public InputError { public: using std::runtime_error::runtime_error; + virtual ~ParseError() throw(); const char * epType() const { return "ParseError"; } virtual const char * epWhat() const { return what(); } }; @@ -45,6 +55,7 @@ namespace ExprParser { class DomainError : public std::domain_error, public InputError { public: using std::domain_error::domain_error; + virtual ~DomainError() throw(); const char * epType() const { return "DomainError"; } virtual const char * epWhat() const { return what(); } }; @@ -55,6 +66,7 @@ namespace ExprParser { class LogicError : public std::logic_error, public Exception { public: using std::logic_error::logic_error; + virtual ~LogicError() throw(); const char * epType() const { return "ImplementationError"; } virtual const char * epWhat() const { return what(); } }; diff --git a/src/simplebuild_dgcode/data/pkgs/Utils/ExprParser/libsrc/Exception.cc b/src/simplebuild_dgcode/data/pkgs/Utils/ExprParser/libsrc/Exception.cc index 811ee327..7a619718 100644 --- a/src/simplebuild_dgcode/data/pkgs/Utils/ExprParser/libsrc/Exception.cc +++ b/src/simplebuild_dgcode/data/pkgs/Utils/ExprParser/libsrc/Exception.cc @@ -1 +1,9 @@ #include "ExprParser/Exception.hh" + +//Out-of-line destructors, so the type information of the exception classes is +//only present in this library (see Exception.hh): +ExprParser::Exception::~Exception() throw() {} +ExprParser::InputError::~InputError() throw() {} +ExprParser::ParseError::~ParseError() throw() {} +ExprParser::DomainError::~DomainError() throw() {} +ExprParser::LogicError::~LogicError() throw() {} diff --git a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/app_mkfile/main.cc b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/app_mkfile/main.cc new file mode 100644 index 00000000..f6793af4 --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/app_mkfile/main.cc @@ -0,0 +1,40 @@ +//Create MCPL files for the tests of G4MCPLPlugins (the stat:sum handling in +//MCPLGen and MCPLWriter). Usage: +// +// sb_g4mcpltests_mkfile FILENAME NPARTICLES WEIGHT INITIAL_RAY_COUNT +// +//A value of "none" for INITIAL_RAY_COUNT means that the file does not get a +//stat:sum:initial_ray_count, while "notfinal" leaves the value at -1 (as when +//the application writing the file was interrupted). + +#include "mcpl.h" +#include +#include +#include + +int main( int argc, char** argv ) +{ + if ( argc != 5 ) { + printf("Usage: %s FILENAME NPARTICLES WEIGHT INITIAL_RAY_COUNT\n",argv[0]); + return 1; + } + const char * irc = argv[4]; + mcpl_outfile_t f = mcpl_create_outfile( argv[1] ); + mcpl_hdr_set_srcname( f, "sb_g4mcpltests_mkfile" ); + if ( std::strcmp( irc, "none" ) != 0 ) + mcpl_hdr_add_stat_sum( f, "initial_ray_count", + std::strcmp( irc, "notfinal" ) == 0 ? -1.0 : std::atof( irc ) ); + mcpl_particle_t * p = mcpl_get_empty_particle( f ); + const int n = std::atoi( argv[2] ); + for ( int i = 0; i < n; ++i ) { + p->pdgcode = 2112; + p->ekin = 0.025e-6 * ( 1 + i );//MeV + p->direction[0] = p->direction[1] = 0.0; + p->direction[2] = 1.0; + p->position[0] = 0.1 * i; + p->weight = std::atof( argv[3] ); + mcpl_add_particle( f, p ); + } + mcpl_close_outfile( f ); + return 0; +} diff --git a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/pkg.info b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/pkg.info index 2d7b0653..9971f728 100644 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/pkg.info +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/pkg.info @@ -1,4 +1,4 @@ -package(USEPKG G4Launcher G4MCPLPlugins MCPLTestData G4StdGeometries G4StdGenerators USEEXT MCPL) +package(USEPKG G4Launcher G4MCPLPlugins MCPLExtra MCPLTestData G4StdGeometries G4StdGenerators USEEXT MCPL) ########################################################## diff --git a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/mcplstatsums b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/mcplstatsums new file mode 100755 index 00000000..4a28af97 --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/mcplstatsums @@ -0,0 +1,24 @@ +#!/usr/bin/env python3 + +#Simulation used by testmcplstatsums: the particles of an MCPL file in an empty +#world (use --mcpl="grabsrc to FILE" to capture them). + +import G4StdGeometries.GeoEmptyWorld as Geo +geo = Geo.create() +geo.material = "Vacuum" +geo.dimension_cm = 200 + +import G4MCPLPlugins.MCPLGen as Gen +gen = Gen.create() + +#A custom stat:sum value for the output files (if requested by the test): +import os +if os.environ.get('MCPLSTATSUMS_TEST_USERVALUE'): + import G4MCPLPlugins.MCPLWriter as Writer + Writer.addStatSum( 'my_beam_time', float(os.environ['MCPLSTATSUMS_TEST_USERVALUE']) ) + +import G4Launcher +g4 = G4Launcher(geo,gen) +g4.setPhysicsList("PL_Empty") +g4.setOutput("none") +g4.go() diff --git a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplgen.log b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplgen.log index 66a47a37..f7464f26 100644 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplgen.log +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplgen.log @@ -29,6 +29,7 @@ G4Launcher:: --> [dbl] roty_degree = 90 G4Launcher:: --> [dbl] rotz_degree = 0 G4Launcher:: --> [int] skip_events = 0 G4Launcher:: --> [flg] allow_zero_weight = no +G4Launcher:: --> [flg] normalise_by_initial_ray_count = yes G4Launcher:: Requested to not set output of GRIFF files: G4Launcher:: Pre-init done G4Launcher:: Calling G4RunManager::Initialize() diff --git a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums new file mode 100755 index 00000000..d6bc2ecb --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums @@ -0,0 +1,126 @@ +#!/usr/bin/env python3 + +#Tests of the stat:sum handling of MCPLGen and MCPLWriter: +# +# * MCPLGen divides the weights by stat:sum:initial_ray_count of the input +# file (if present), unless normalise_by_initial_ray_count is disabled. +# * MCPLWriter writes stat:sum:geant4_nevts and stat:sum:initial_geant4_weights +# (the number of simulated events and the sum of the weights of their +# primary particles), which are added up when the output of multiple +# processes is merged. +# * MCPLWriter records stat:sum:initial_ray_count of the input file as +# stat:sum:original_initial_ray_count, with the weights as used in Geant4. +# With the opt_irc keyword, it is instead recorded as initial_ray_count, with +# the weights in the same convention as in the input file. +# * Values describing the whole job (like the above) are divided between the +# processes when multiprocessing, so the merged file gets the value. +# * Custom stat:sum values can be added (here from Python). +# * sb_mcplextra_filterfile marks stat:sum values as unavailable when only a +# part of the input file is used. + +import subprocess +import mcpl + +def mkfile( fn, nparticles, weight, initial_ray_count ): + subprocess.run( [ 'sb_g4mcpltests_mkfile', fn, str(nparticles), str(weight), + str(initial_ray_count) ], check = True ) + +def sim( title, infile, *args, expect_ok = True, env_extra = None, mcplkw = 'grabsrc' ): + print() + print( f'==> {title}' ) + outfile = 'out.mcpl' + cmd = [ 'sb_g4mcpltests_mcplstatsums', f'input_file={infile}', + f'--mcpl={mcplkw} to {outfile}' ] + list(args) + import os + env = dict( os.environ ) + env.update( env_extra or {} ) + r = subprocess.run( cmd, capture_output = True, text = True, check = False, env = env ) + #Only show the messages of MCPLGen and MCPLWriter (others are not relevant + #here). With multiple processes, a message might start in the middle of a + #line of output from another process, so look for them anywhere in a line: + for line in ( r.stdout + r.stderr ).splitlines(): + for m in ( 'MCPLGen: Dividing', 'MCPLGen: Using', 'MCPLGen ERROR', + 'MCPLWriter: Done merging', 'MCPLWriter: Writing stat:sum' ): + if m in line: + print( ' ', line[line.index(m):] ) + break + if not expect_ok: + assert r.returncode != 0, 'simulation unexpectedly succeeded' + print( ' (simulation failed as expected)' ) + return + assert r.returncode == 0, 'simulation failed:\n' + r.stdout + r.stderr + f = mcpl.MCPLFile( outfile + '.gz' ) + weights = sorted( float(p.weight) for p in f.particles ) + print( f' Output weights: {weights}' ) + for k, v in sorted( f.stat_sum.items() ): + print( f' Output stat:sum:{k} = {"" if v is None else format(v,"g")}' ) + +mkfile( 'rays.mcpl', 5, 1000.0, 1000 ) +mkfile( 'norays.mcpl', 5, 2.0, 'none' ) +mkfile( 'notfinal.mcpl', 5, 1000.0, 'notfinal' ) +mkfile( 'zero.mcpl', 5, 1000.0, 0 ) + +sim( 'File with initial_ray_count', 'rays.mcpl' ) +sim( 'File with initial_ray_count, not normalised', 'rays.mcpl', + 'normalise_by_initial_ray_count=0' ) +sim( 'File with initial_ray_count, skipping events', 'rays.mcpl', 'skip_events=2' ) +sim( 'File with initial_ray_count, limited number of events', 'rays.mcpl', '-n2' ) +sim( 'File with initial_ray_count, multiple processes', 'rays.mcpl', '-j2' ) +sim( 'File with initial_ray_count, multiple processes, limited number of events', + 'rays.mcpl', '-j2', '-n3' ) +sim( 'File without initial_ray_count', 'norays.mcpl' ) +sim( 'File with initial_ray_count which is not finalised', 'notfinal.mcpl', + expect_ok = False ) +sim( 'File with initial_ray_count which is not finalised, not normalised', + 'notfinal.mcpl', 'normalise_by_initial_ray_count=0' ) +sim( 'File with initial_ray_count of zero', 'zero.mcpl', expect_ok = False ) + +#Recording the initial_ray_count of the input file as initial_ray_count, with +#the weights in the same convention as in the input file (opt_irc): +sim( 'File with initial_ray_count, opt_irc', 'rays.mcpl', mcplkw = 'grabsrc,opt_irc' ) +sim( 'File with initial_ray_count, opt_irc, not normalised', 'rays.mcpl', + 'normalise_by_initial_ray_count=0', mcplkw = 'grabsrc,opt_irc' ) +sim( 'File with initial_ray_count, opt_irc, multiple processes', 'rays.mcpl', '-j2', + mcplkw = 'grabsrc,opt_irc' ) +sim( 'File without initial_ray_count, opt_irc', 'norays.mcpl', mcplkw = 'grabsrc,opt_irc', + expect_ok = False ) + +#Reading the output files again. With opt_irc, the weights are normalised +#again, so they are the same in Geant4 as in the original simulation: +sim( 'File with initial_ray_count, opt_irc, for reading again', 'rays.mcpl', + mcplkw = 'grabsrc,opt_irc' ) +subprocess.run( [ 'mv', 'out.mcpl.gz', 'again.mcpl.gz' ], check = True ) +sim( 'Reading the output again', 'again.mcpl.gz' ) +#By default, the output has the weights as used in Geant4, and no +#initial_ray_count (so it is not normalised again): +subprocess.run( [ 'mv', 'out.mcpl.gz', 'again2.mcpl.gz' ], check = True ) +sim( 'Reading the output of that again', 'again2.mcpl.gz' ) + +#Custom stat:sum values added by user code (describing the whole simulation, so +#they are only added once when merging the output of multiple processes): +sim( 'Custom stat:sum value', 'rays.mcpl', env_extra = { 'MCPLSTATSUMS_TEST_USERVALUE' : '17.5' } ) +sim( 'Custom stat:sum value, multiple processes', 'rays.mcpl', '-j2', + env_extra = { 'MCPLSTATSUMS_TEST_USERVALUE' : '17.5' } ) + +#sb_mcplextra_filterfile keeps the stat:sum values when filtering the particles, +#but marks them as unavailable when only a part of the file is used due to a +#limit on the number of particles: +def filterfile( title, *args ): + print() + print( f'==> {title}' ) + r = subprocess.run( [ 'sb_mcplextra_filterfile', '-n', 'rays.mcpl', 'filtered.mcpl' ] + list(args), + capture_output = True, text = True, check = False ) + for line in r.stdout.splitlines(): + if 'stat:sum' in line: + print( ' ', line ) + assert r.returncode == 0, r.stdout + r.stderr + f = mcpl.MCPLFile( 'filtered.mcpl' ) + print( f' Output particles: {f.nparticles}' ) + for k, v in sorted( f.stat_sum.items() ): + print( f' Output stat:sum:{k} = {"" if v is None else format(v,"g")}' ) + +filterfile( 'Filtering a file', 'x > 0.15cm' ) +filterfile( 'Filtering a file with a limit which is not reached', '-l10' ) +filterfile( 'Filtering a file with a limit', '-l2' ) +filterfile( 'Filtering a file with a limit, with a filter', '-l2', 'x > 0.15cm' ) +sim( 'File where initial_ray_count is unavailable', 'filtered.mcpl', expect_ok = False ) diff --git a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log new file mode 100644 index 00000000..f6e001e3 --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log @@ -0,0 +1,148 @@ + +==> File with initial_ray_count + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + Output weights: [1.0, 1.0, 1.0, 1.0, 1.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5 + Output stat:sum:original_initial_ray_count = 1000 + +==> File with initial_ray_count, not normalised + MCPLGen: Using the weights in the file as stored (normalise_by_initial_ray_count is disabled, stat:sum:initial_ray_count is 1000) + Output weights: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5000 + Output stat:sum:original_initial_ray_count = 1000 + +==> File with initial_ray_count, skipping events + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + Output weights: [1.0, 1.0, 1.0] + Output stat:sum:geant4_nevts = 3 + Output stat:sum:initial_geant4_weights = 3 + Output stat:sum:original_initial_ray_count = 1000 + +==> File with initial_ray_count, limited number of events + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + Output weights: [1.0, 1.0] + Output stat:sum:geant4_nevts = 2 + Output stat:sum:initial_geant4_weights = 2 + Output stat:sum:original_initial_ray_count = 1000 + +==> File with initial_ray_count, multiple processes + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + MCPLWriter: Done merging output from 2 processes into file out.mcpl + Output weights: [1.0, 1.0, 1.0, 1.0, 1.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5 + Output stat:sum:original_initial_ray_count = 1000 + +==> File with initial_ray_count, multiple processes, limited number of events + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + MCPLWriter: Done merging output from 2 processes into file out.mcpl + Output weights: [1.0, 1.0, 1.0] + Output stat:sum:geant4_nevts = 3 + Output stat:sum:initial_geant4_weights = 3 + Output stat:sum:original_initial_ray_count = 1000 + +==> File without initial_ray_count + Output weights: [2.0, 2.0, 2.0, 2.0, 2.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 10 + +==> File with initial_ray_count which is not finalised + MCPLGen ERROR - The input file has a stat:sum:initial_ray_count value which is not available (for instance due to an interrupted run of the application writing the file, or because the file contains only a part of the particles of another file), so the weights can not be normalised. Set normalise_by_initial_ray_count=false to use the weights as stored in the file. + (simulation failed as expected) + +==> File with initial_ray_count which is not finalised, not normalised + MCPLGen: Using the weights in the file as stored (normalise_by_initial_ray_count is disabled, stat:sum:initial_ray_count is -1) + Output weights: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5000 + Output stat:sum:original_initial_ray_count = + +==> File with initial_ray_count of zero + MCPLGen ERROR - The input file has an invalid stat:sum:initial_ray_count value (0), so the weights can not be normalised. Set normalise_by_initial_ray_count=false to use the weights as stored in the file. + (simulation failed as expected) + +==> File with initial_ray_count, opt_irc + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + Output weights: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5000 + Output stat:sum:initial_ray_count = 1000 + +==> File with initial_ray_count, opt_irc, not normalised + MCPLGen: Using the weights in the file as stored (normalise_by_initial_ray_count is disabled, stat:sum:initial_ray_count is 1000) + Output weights: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5000 + Output stat:sum:initial_ray_count = 1000 + +==> File with initial_ray_count, opt_irc, multiple processes + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + MCPLWriter: Done merging output from 2 processes into file out.mcpl + Output weights: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5000 + Output stat:sum:initial_ray_count = 1000 + +==> File without initial_ray_count, opt_irc + MCPLWriter: Writing stat:sum:initial_ray_count was requested, but there is no MCPL input file (read with MCPLGen) with a valid stat:sum:initial_ray_count value + (simulation failed as expected) + +==> File with initial_ray_count, opt_irc, for reading again + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + Output weights: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5000 + Output stat:sum:initial_ray_count = 1000 + +==> Reading the output again + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + Output weights: [1.0, 1.0, 1.0, 1.0, 1.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5 + Output stat:sum:original_initial_ray_count = 1000 + +==> Reading the output of that again + Output weights: [1.0, 1.0, 1.0, 1.0, 1.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5 + +==> Custom stat:sum value + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + Output weights: [1.0, 1.0, 1.0, 1.0, 1.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5 + Output stat:sum:my_beam_time = 17.5 + Output stat:sum:original_initial_ray_count = 1000 + +==> Custom stat:sum value, multiple processes + MCPLGen: Dividing all weights by stat:sum:initial_ray_count = 1000 (disable with normalise_by_initial_ray_count=false) + MCPLWriter: Done merging output from 2 processes into file out.mcpl + Output weights: [1.0, 1.0, 1.0, 1.0, 1.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5 + Output stat:sum:my_beam_time = 17.5 + Output stat:sum:original_initial_ray_count = 1000 + +==> Filtering a file + Output particles: 3 + Output stat:sum:initial_ray_count = 1000 + +==> Filtering a file with a limit which is not reached + Output particles: 5 + Output stat:sum:initial_ray_count = 1000 + +==> Filtering a file with a limit + Marking stat:sum values as unavailable, since the limit on the number of particles means that only a part of the input file was used. + Output particles: 2 + Output stat:sum:initial_ray_count = + +==> Filtering a file with a limit, with a filter + Marking stat:sum values as unavailable, since the limit on the number of particles means that only a part of the input file was used. + Output particles: 2 + Output stat:sum:initial_ray_count = + +==> File where initial_ray_count is unavailable + MCPLGen ERROR - The input file has a stat:sum:initial_ray_count value which is not available (for instance due to an interrupted run of the application writing the file, or because the file contains only a part of the particles of another file), so the weights can not be normalised. Set normalise_by_initial_ray_count=false to use the weights as stored in the file. + (simulation failed as expected) diff --git a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplwrite.log b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplwrite.log index d4e7fc5b..51f8ebce 100644 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplwrite.log +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplwrite.log @@ -871,12 +871,12 @@ Opened MCPL file target.mcpl.gz: Basic info Format : MCPL-3 No. of particles : 0 - Header storage : 931 bytes + Header storage : 1325 bytes Data storage : 0 bytes Custom meta data Source : "Geant4" - Number of comments : 8 + Number of comments : 11 -> comment 0 : "Created with the Geant4 MCPLWriter in the ESS/dgcode framework" -> comment 1 : "MPCLWriter volumes considered : ['Target']" -> comment 2 : "MPCLWriter steps considered : " @@ -885,6 +885,9 @@ Opened MCPL file target.mcpl.gz: -> comment 5 : "MPCLWriter track kill strategy : " -> comment 6 : "ESS/dgcode geometry module : G4StdGeometries/GeoSlab" -> comment 7 : "ESS/dgcode generator module : G4StdGenerators/SimpleGen" + -> comment 8 : "MCPLWriter stat:sum values : geant4_nevts is the number of simulated events, initial_geant4_weights is the sum of the weights of their primary particles (in the convention of the weights in this file), original_initial_ray_count is the initial_ray_count of the input file (if any)" + -> comment 9 : "stat:sum:geant4_nevts: 10" + -> comment 10 : "stat:sum:initial_geant4_weights: 10" Number of blobs : 2 -> 106 bytes of data with key "ESS/dgcode_geopars" -> 231 bytes of data with key "ESS/dgcode_genpars" @@ -903,12 +906,12 @@ Opened MCPL file recordfwd.mcpl.gz: Basic info Format : MCPL-3 No. of particles : 6 - Header storage : 858 bytes + Header storage : 1252 bytes Data storage : 192 bytes Custom meta data Source : "Geant4" - Number of comments : 8 + Number of comments : 11 -> comment 0 : "Created with the Geant4 MCPLWriter in the ESS/dgcode framework" -> comment 1 : "MPCLWriter volumes considered : ['RecordFwd']" -> comment 2 : "MPCLWriter steps considered : " @@ -917,6 +920,9 @@ Opened MCPL file recordfwd.mcpl.gz: -> comment 5 : "MPCLWriter track kill strategy : " -> comment 6 : "ESS/dgcode geometry module : G4StdGeometries/GeoSlab" -> comment 7 : "ESS/dgcode generator module : G4StdGenerators/SimpleGen" + -> comment 8 : "MCPLWriter stat:sum values : geant4_nevts is the number of simulated events, initial_geant4_weights is the sum of the weights of their primary particles (in the convention of the weights in this file), original_initial_ray_count is the initial_ray_count of the input file (if any)" + -> comment 9 : "stat:sum:geant4_nevts: 10" + -> comment 10 : "stat:sum:initial_geant4_weights: 10" Number of blobs : 2 -> 106 bytes of data with key "ESS/dgcode_geopars" -> 231 bytes of data with key "ESS/dgcode_genpars"