From 6a23688eb0472d928eda49445255b1940baf1313 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Thu, 5 Nov 2026 22:51:46 +0000 Subject: [PATCH 1/8] MCPL: normalise MCPLGen weights by initial_ray_count, add stat sums to MCPLWriter output - MCPLGen divides the weights by stat:sum:initial_ray_count if the input file has it (as files from the McStas MCPL_output component with weight_mode=1, the default since McStas 3.10, do). The value of the whole file is always used, so results of jobs processing parts of a file can be added up. Files without the value are used as before, and files with a value which is not finalised or not positive give an error. The new parameter normalise_by_initial_ray_count (default true) can be disabled to always use the weights as stored. - MCPLWriter writes stat:sum:geant4_nevts (the number of simulated events) and stat:sum:initial_geant4_weights (the sum of the weights of their primary particles) to the output files. These are collected via a post-generation callback on the generator, and are added up when files are merged. - Require MCPL 2.2.8 or later. Assisted-by: Claude --- .github/resources/conda-setup.yml | 2 +- doc/source/mcpl.rst | 25 ++++++ pyproject.toml | 4 +- .../MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc | 59 +++++++++++++- .../G4MCPLPlugins/pymod__MCPLWriter/mod.cc | 81 ++++++++++++++++++- .../UnitTests/G4MCPLTests/app_mkfile/main.cc | 40 +++++++++ .../G4MCPLTests/scripts/mcplstatsums | 18 +++++ .../G4MCPLTests/scripts/testmcplgen.log | 1 + .../G4MCPLTests/scripts/testmcplstatsums | 63 +++++++++++++++ .../G4MCPLTests/scripts/testmcplstatsums.log | 61 ++++++++++++++ .../G4MCPLTests/scripts/testmcplwrite.log | 14 +++- 11 files changed, 353 insertions(+), 15 deletions(-) create mode 100644 src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/app_mkfile/main.cc create mode 100755 src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/mcplstatsums create mode 100755 src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums create mode 100644 src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log diff --git a/.github/resources/conda-setup.yml b/.github/resources/conda-setup.yml index c2d6a38..16cbac7 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/doc/source/mcpl.rst b/doc/source/mcpl.rst index 1ebe927..c1f31ea 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,19 @@ options, by using the module. See an example :sbpkg:`here ` for how to use it. +The MCPL files written by Geant4 contain two ``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. + +The weights of the particles in the file are the weights in the Geant4 +simulation (i.e. any ``stat:sum:initial_ray_count`` of an input file read with +MCPLGen is already accounted for, and is not copied to the output). + Modify MCPL files ----------------- diff --git a/pyproject.toml b/pyproject.toml index f1db4b6..a36c8ad 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/MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc index 3bd7ea9..8356eae 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 @@ -4,6 +4,9 @@ #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 +15,7 @@ #include "CLHEP/Geometry/Vector3D.h" #include "CLHEP/Geometry/Transform3D.h" #include +#include #include #include #include "G4ParticleTable.hh" @@ -33,6 +37,7 @@ class MCPLGen : public G4Interfaces::ParticleGenBase { void skip_forward(unsigned nskip); void delayed_init(); void setPDG(int); + void init_weight_scale(); mcpl_file_t m_mcplfile; const mcpl_particle_t * m_p;//particle to be used for next generation bool m_unfiltered; @@ -48,6 +53,7 @@ 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 }; MCPLGen::MCPLGen() @@ -58,7 +64,8 @@ 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_mcplfile.internal = 0; @@ -83,6 +90,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() @@ -171,6 +186,8 @@ void MCPLGen::delayed_init() m_has_polarisation = mcpl_hdr_has_polarisation(m_mcplfile); + init_weight_scale(); + unsigned nskip = getParameterInt("skip_events"); if (nskip) skip_forward(nskip); @@ -191,6 +208,38 @@ void MCPLGen::delayed_init() } } +void MCPLGen::init_weight_scale() +{ + m_weight_scale = 1.0; + const double nrays = mcpl_hdr_stat_sum(m_mcplfile,"initial_ray_count"); + //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); + 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 was not" + " finalised (for instance due to an interrupted run of the application writing the 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 = 1.0 / nrays; + if (doprint) + printf("MCPLGen: Dividing all weights by stat:sum:initial_ray_count = %g" + " (disable with normalise_by_initial_ray_count=false)\n",nrays); +} + void MCPLGen::setPDG(int p) { assert(m_lastpdgcode!=p); @@ -273,11 +322,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 2b813b0..ed46197 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 @@ -7,8 +7,14 @@ #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" @@ -34,6 +40,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) { @@ -73,7 +122,8 @@ namespace G4MCPLWriter { 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 +139,8 @@ 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) { std::memset(&m_p,0,sizeof(m_p)); } @@ -125,6 +176,11 @@ namespace G4MCPLWriter { return; m_closed = true; assert(m_f.internal); + if (m_evtstats) { + //Final values of the stat:sum entries reserved in initmcpl(): + mcpl_hdr_add_stat_sum(m_f,"geant4_nevts",(double)m_evtstats->nevts()); + mcpl_hdr_add_stat_sum(m_f,"initial_geant4_weights",m_evtstats->sumWeights()); + } if (!FrameworkGlobals::isForked()) mcpl_closeandgzip_outfile(m_f); else @@ -268,6 +324,13 @@ 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); + } + 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); @@ -290,6 +353,7 @@ namespace G4MCPLWriter { mcpl_particle_t m_p; MCPLOutputMerger * m_om; std::vector> m_comments_and_blobs; + EventStatsPtr m_evtstats; }; void MCPLOutputMerger::merge() { @@ -370,6 +434,10 @@ 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!"); @@ -448,7 +516,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 +605,9 @@ 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"); if (geo) { char * geodata(0); @@ -588,7 +660,8 @@ namespace G4MCPLWriter { m_expr_builder, m_eval_filter, m_eval_flags, - createCommentsAndBlobs()); + createCommentsAndBlobs(), + m_evtstats); G4SDManager::GetSDMpointer()->AddNewDetector(m_sd); } auto currentsd = lv->GetSensitiveDetector(); 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 0000000..f6793af --- /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/scripts/mcplstatsums b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/mcplstatsums new file mode 100755 index 0000000..9bc9e74 --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/mcplstatsums @@ -0,0 +1,18 @@ +#!/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() + +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 66a47a3..f7464f2 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 0000000..bb51fad --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums @@ -0,0 +1,63 @@ +#!/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. + +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 ): + print() + print( f'==> {title}' ) + outfile = 'out.mcpl' + cmd = [ 'sb_g4mcpltests_mcplstatsums', f'input_file={infile}', + f'--mcpl=grabsrc to {outfile}' ] + list(args) + r = subprocess.run( cmd, capture_output = True, text = True, check = False ) + #Only show the lines of MCPLGen and MCPLWriter (others are not relevant here): + for line in ( r.stdout + r.stderr ).splitlines(): + if line.startswith( ( 'MCPLGen', 'MCPLWriter: Done merging' ) ): + print( ' ', line ) + 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} = {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 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 ) + +#Reading a file written by MCPLWriter again does not normalise its weights a +#second time (since it has no initial_ray_count): +sim( 'File with initial_ray_count, for reading again', 'rays.mcpl' ) +subprocess.run( [ 'mv', 'out.mcpl.gz', 'again.mcpl.gz' ], check = True ) +sim( 'Reading the output again', 'again.mcpl.gz' ) 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 0000000..8e64a3e --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log @@ -0,0 +1,61 @@ + +==> 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 + +==> 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 + +==> 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 + +==> 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 + +==> 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 + +==> 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 was not finalised (for instance due to an interrupted run of the application writing the 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 + +==> 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, for reading 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 + +==> Reading the output 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 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 d4e7fc5..774cf30 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 : 1197 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" + -> 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 : 1124 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" + -> 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" From 1f7453fe31848f36256ec1c0f91022e8c9dfefb2 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Fri, 6 Nov 2026 06:39:32 +0000 Subject: [PATCH 2/8] MCPLExtra: mark stat:sum values unavailable when filterfile truncates a file sb_mcplextra_filterfile keeps stat:sum values when filtering particles, but when only a part of the input file is used due to the -l limit, the values are marked as unavailable (the file was split rather than filtered). MCPLGen refuses to normalise by such values by default. Assisted-by: Claude --- doc/source/mcpl.rst | 6 +++++ .../MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc | 9 +++---- .../MCPL/MCPLExtra/app_filterfile/main.cc | 18 +++++++++++++- .../pkgs_val/UnitTests/G4MCPLTests/pkg.info | 2 +- .../G4MCPLTests/scripts/testmcplstatsums | 23 ++++++++++++++++++ .../G4MCPLTests/scripts/testmcplstatsums.log | 24 ++++++++++++++++++- 6 files changed, 75 insertions(+), 7 deletions(-) diff --git a/doc/source/mcpl.rst b/doc/source/mcpl.rst index c1f31ea..5e8e601 100644 --- a/doc/source/mcpl.rst +++ b/doc/source/mcpl.rst @@ -139,6 +139,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/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc index 8356eae..1ed7314 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 @@ -224,10 +224,11 @@ void MCPLGen::init_weight_scale() 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 was not" - " finalised (for instance due to an interrupted run of the application writing the 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"); + 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" 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 12405d9..4d94fd9 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_val/UnitTests/G4MCPLTests/pkg.info b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/pkg.info index 2d7b065..9971f72 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/testmcplstatsums b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums index bb51fad..5488709 100755 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums @@ -61,3 +61,26 @@ sim( 'File with initial_ray_count of zero', 'zero.mcpl', expect_ok = False ) sim( 'File with initial_ray_count, for reading again', 'rays.mcpl' ) subprocess.run( [ 'mv', 'out.mcpl.gz', 'again.mcpl.gz' ], check = True ) sim( 'Reading the output again', 'again.mcpl.gz' ) + +#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} = {v}' ) + +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 index 8e64a3e..578a37b 100644 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log @@ -36,7 +36,7 @@ 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 was not finalised (for instance due to an interrupted run of the application writing the file), so the weights can not be normalised. Set normalise_by_initial_ray_count=false to use the weights as stored in the file. + 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 @@ -59,3 +59,25 @@ 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 + +==> Filtering a file + Output particles: 3 + Output stat:sum:initial_ray_count = 1000.0 + +==> Filtering a file with a limit which is not reached + Output particles: 5 + Output stat:sum:initial_ray_count = 1000.0 + +==> 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 = None + +==> 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 = None + +==> 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) From dc24b3613db4030b3b89c58742c80da24f352940 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Fri, 6 Nov 2026 06:54:25 +0000 Subject: [PATCH 3/8] MCPLWriter: transfer stat:sum:initial_ray_count of MCPLGen input files If the input file of MCPLGen has a stat:sum:initial_ray_count value, it is transferred to the output files of MCPLWriter, with the weights written in the same convention as in the input file (i.e. multiplied by the value when MCPLGen divided the weights by it), so reading the output again gives the same weights in Geant4. The value is marked as unavailable unless all particles of the input file were used. When multiprocessing, the value is only set in the file of the parent process (the files of the child processes get 0), since the values are added up when the files are merged. The information is provided by MCPLGen via the new G4MCPL/InitialRayCountInfo.hh interface. Assisted-by: Claude --- doc/source/mcpl.rst | 10 +++- .../MCPL/G4MCPL/libinc/InitialRayCountInfo.hh | 36 ++++++++++++ .../MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc | 57 ++++++++++++++++--- .../G4MCPLPlugins/pymod__MCPLWriter/mod.cc | 46 +++++++++++++-- .../G4MCPLTests/scripts/testmcplstatsums | 17 ++++-- .../G4MCPLTests/scripts/testmcplstatsums.log | 37 ++++++++---- 6 files changed, 174 insertions(+), 29 deletions(-) create mode 100644 src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/InitialRayCountInfo.hh diff --git a/doc/source/mcpl.rst b/doc/source/mcpl.rst index 5e8e601..e0ceb29 100644 --- a/doc/source/mcpl.rst +++ b/doc/source/mcpl.rst @@ -115,9 +115,13 @@ or when merging the output of multiple processes): * ``stat:sum:initial_geant4_weights``: The sum of the weights of the primary particles of the simulated events. -The weights of the particles in the file are the weights in the Geant4 -simulation (i.e. any ``stat:sum:initial_ray_count`` of an input file read with -MCPLGen is already accounted for, and is not copied to the output). +If the input file read with MCPLGen has a ``stat:sum:initial_ray_count`` value, +it is transferred to the output file, and the weights are written in the same +convention as in the input file (i.e. multiplied by the value if MCPLGen divided +the weights by it). So reading the output file with MCPLGen again gives the same +weights as in the original simulation. Since the value is only valid when all +particles of the input file were used, it is marked as unavailable when this +was not the case (e.g. with ``skip_events`` or a limited number of events). Modify MCPL files ----------------- diff --git a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/InitialRayCountInfo.hh b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/InitialRayCountInfo.hh new file mode 100644 index 0000000..17ce857 --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/InitialRayCountInfo.hh @@ -0,0 +1,36 @@ +#ifndef G4MCPL_InitialRayCountInfo_hh +#define G4MCPL_InitialRayCountInfo_hh + +namespace G4MCPL { + + //Interface implemented by particle generators reading particles from MCPL + //files (i.e. G4MCPLPlugins/MCPLGen), with information about the + //stat:sum:initial_ray_count value of the input file. This is used by + //G4MCPLPlugins/MCPLWriter to transfer the value to its output files, with + //the weights written in the same convention as in the input file (i.e. the + //weights of the input file divided by the value are the weights per initial + //ray). + + class InitialRayCountInfo { + public: + virtual ~InitialRayCountInfo() = default; + + //The stat:sum:initial_ray_count value of the input file: -2 if not + //present, -1 if not available, and otherwise the (positive) value. This + //only reads the header of the input file, so it can be called at any time + //after the generator is configured: + virtual double inputInitialRayCount() const = 0; + + //The factor applied to the weights of the input file when generating + //particles (1/initial_ray_count when normalising the weights, otherwise + //1): + virtual double inputWeightFactor() const = 0; + + //Whether all particles of the input file (in this process, when + //multiprocessing) were used: + virtual bool usedAllInputParticles() const = 0; + }; + +} + +#endif 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 1ed7314..af46f6e 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,4 +1,5 @@ #include "G4MCPL/G4MCPLUserFlags.hh" +#include "G4MCPL/InitialRayCountInfo.hh" #include "G4Interfaces/ParticleGenToPy.hh" #include "G4Interfaces/ParticleGenBase.hh" #include "G4Interfaces/FrameworkGlobals.hh" @@ -21,13 +22,19 @@ #include "G4ParticleTable.hh" #include "G4IonTable.hh" -class MCPLGen : public G4Interfaces::ParticleGenBase { +class MCPLGen : public G4Interfaces::ParticleGenBase, + public G4MCPL::InitialRayCountInfo { public: MCPLGen(); virtual ~MCPLGen(); void init(); void gen(G4Event*); + //For MCPLWriter (see G4MCPL/InitialRayCountInfo.hh): + double inputInitialRayCount() const override; + double inputWeightFactor() const override; + bool usedAllInputParticles() const override; + virtual bool unlimited() const { return false; }//only has the events available in the input protected: @@ -38,6 +45,7 @@ class MCPLGen : public G4Interfaces::ParticleGenBase { void delayed_init(); void setPDG(int); void init_weight_scale(); + std::string resolved_input_file() const; mcpl_file_t m_mcplfile; const mcpl_particle_t * m_p;//particle to be used for next generation bool m_unfiltered; @@ -54,6 +62,7 @@ class MCPLGen : public G4Interfaces::ParticleGenBase { 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) }; MCPLGen::MCPLGen() @@ -65,7 +74,8 @@ MCPLGen::MCPLGen() m_lastpdgcode(0), m_has_polarisation(false), m_allow_zero_weight(false), - m_weight_scale(1.0) + m_weight_scale(1.0), + m_input_irc(std::nan("")) { m_mcplfile.internal = 0; @@ -172,17 +182,50 @@ 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; +} + +bool MCPLGen::usedAllInputParticles() const +{ + return getParameterInt("skip_events") == 0 && reachedLimit(); +} + +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); @@ -211,7 +254,7 @@ void MCPLGen::delayed_init() void MCPLGen::init_weight_scale() { m_weight_scale = 1.0; - const double nrays = mcpl_hdr_stat_sum(m_mcplfile,"initial_ray_count"); + const double nrays = inputInitialRayCount(); //Only print in one process when multiprocessing: const bool doprint = FrameworkGlobals::isParent(); if (!getParameterBoolean("normalise_by_initial_ray_count")) { @@ -235,7 +278,7 @@ void MCPLGen::init_weight_scale() " 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 = 1.0 / nrays; + 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); 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 ed46197..3bc1221 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,6 +2,7 @@ #include "G4Interfaces/ComponentHandles.hh" #include "G4MCPL/G4MCPLUserFlags.hh" +#include "G4MCPL/InitialRayCountInfo.hh" #include "G4Interfaces/FrameworkGlobals.hh" #include "G4Interfaces/GeoConstructBase.hh" #include "G4Interfaces/ParticleGenBase.hh" @@ -123,7 +124,8 @@ namespace G4MCPLWriter { EP::Evaluator filter_evaluator, EP::Evaluator flag_evaluator, std::vector> comments_and_blobs, - EventStatsPtr evtstats) + EventStatsPtr evtstats, + const G4MCPL::InitialRayCountInfo* ircinfo) : G4VSensitiveDetector(name), m_expr_builder(builder), m_eval_filter(filter_evaluator), @@ -140,7 +142,10 @@ namespace G4MCPLWriter { m_closed(false), m_om(0), m_comments_and_blobs(comments_and_blobs), - m_evtstats(evtstats) + m_evtstats(evtstats), + m_ircinfo(ircinfo), + m_transfer_irc(false), + m_weight_factor(1.0) { std::memset(&m_p,0,sizeof(m_p)); } @@ -181,6 +186,23 @@ namespace G4MCPLWriter { mcpl_hdr_add_stat_sum(m_f,"geant4_nevts",(double)m_evtstats->nevts()); mcpl_hdr_add_stat_sum(m_f,"initial_geant4_weights",m_evtstats->sumWeights()); } + if (m_transfer_irc) { + //The initial_ray_count of the input file is only valid if all of its + //particles were used. When multiprocessing, the value is only set in + //the file of the parent process (since the values are added up when + //the files are merged), and the files of all processes must have used + //all their particles (-1 in any of the files makes the merged value + //unavailable): + // + //NB: The value would preferably be set only after merging the files, + //but MCPL currently does not support editing stat:sum values of + //existing files. + const double irc = m_ircinfo->inputInitialRayCount(); + double value = -1.0; + if ( irc > 0.0 && m_ircinfo->usedAllInputParticles() ) + value = FrameworkGlobals::isChild() ? 0.0 : irc; + mcpl_hdr_add_stat_sum(m_f,"initial_ray_count",value); + } if (!FrameworkGlobals::isForked()) mcpl_closeandgzip_outfile(m_f); else @@ -251,7 +273,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(); @@ -331,10 +353,20 @@ namespace G4MCPLWriter { mcpl_hdr_add_stat_sum(m_f,"initial_geant4_weights",-1.0); } + //Transfer the stat:sum:initial_ray_count of an input file (read by + //MCPLGen), with the weights in the same convention as in the input file + //(i.e. undoing any normalisation by MCPLGen), so the weights per initial + //ray are the same when reading the file again: + if ( m_ircinfo && m_ircinfo->inputInitialRayCount() != -2.0 ) { + m_transfer_irc = true; + m_weight_factor = 1.0 / m_ircinfo->inputWeightFactor(); + mcpl_hdr_add_stat_sum(m_f,"initial_ray_count",-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() @@ -354,6 +386,9 @@ namespace G4MCPLWriter { MCPLOutputMerger * m_om; std::vector> m_comments_and_blobs; EventStatsPtr m_evtstats; + const G4MCPL::InitialRayCountInfo* m_ircinfo; + bool m_transfer_irc; + double m_weight_factor; }; void MCPLOutputMerger::merge() { @@ -661,7 +696,8 @@ namespace G4MCPLWriter { m_eval_filter, m_eval_flags, createCommentsAndBlobs(), - m_evtstats); + m_evtstats, + dynamic_cast(m_gen)); G4SDManager::GetSDMpointer()->AddNewDetector(m_sd); } auto currentsd = lv->GetSensitiveDetector(); diff --git a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums index 5488709..e2c52c3 100755 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums @@ -8,6 +8,11 @@ # (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 transfers stat:sum:initial_ray_count of the input file to the +# output file (marked as unavailable unless all particles of the input file +# were used), with the weights in the same convention as in the input file. +# * sb_mcplextra_filterfile marks stat:sum values as unavailable when only a +# part of the input file is used. import subprocess import mcpl @@ -36,7 +41,7 @@ def sim( title, infile, *args, expect_ok = True ): 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} = {v:g}' ) + 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' ) @@ -49,6 +54,8 @@ sim( 'File with initial_ray_count, not normalised', 'rays.mcpl', 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 ) @@ -56,8 +63,10 @@ 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 ) -#Reading a file written by MCPLWriter again does not normalise its weights a -#second time (since it has no initial_ray_count): +#The initial_ray_count of the input file is transferred to the output file if +#all particles of the input file were used, with the weights in the same +#convention as in the input file. So reading the output again gives the same +#weights in Geant4: sim( 'File with initial_ray_count, for reading again', 'rays.mcpl' ) subprocess.run( [ 'mv', 'out.mcpl.gz', 'again.mcpl.gz' ], check = True ) sim( 'Reading the output again', 'again.mcpl.gz' ) @@ -77,7 +86,7 @@ def filterfile( title, *args ): 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} = {v}' ) + 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' ) 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 index 578a37b..546ea8a 100644 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log @@ -1,34 +1,47 @@ ==> 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 weights: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] Output stat:sum:geant4_nevts = 5 Output stat:sum:initial_geant4_weights = 5 + Output stat:sum: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: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 weights: [1000.0, 1000.0, 1000.0] Output stat:sum:geant4_nevts = 3 Output stat:sum:initial_geant4_weights = 3 + Output stat:sum:initial_ray_count = ==> 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 weights: [1000.0, 1000.0] Output stat:sum:geant4_nevts = 2 Output stat:sum:initial_geant4_weights = 2 + Output stat:sum:initial_ray_count = ==> 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 weights: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] Output stat:sum:geant4_nevts = 5 Output stat:sum:initial_geant4_weights = 5 + Output stat:sum: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: [1000.0, 1000.0, 1000.0] + Output stat:sum:geant4_nevts = 3 + Output stat:sum:initial_geant4_weights = 3 + Output stat:sum:initial_ray_count = ==> File without initial_ray_count Output weights: [2.0, 2.0, 2.0, 2.0, 2.0] @@ -44,6 +57,7 @@ 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 = ==> 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. @@ -51,32 +65,35 @@ ==> File with initial_ray_count, for reading 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 weights: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] Output stat:sum:geant4_nevts = 5 Output stat:sum:initial_geant4_weights = 5 + Output stat:sum:initial_ray_count = 1000 ==> Reading the output again - Output weights: [1.0, 1.0, 1.0, 1.0, 1.0] + 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 = 5 + Output stat:sum:initial_ray_count = 1000 ==> Filtering a file Output particles: 3 - Output stat:sum:initial_ray_count = 1000.0 + 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.0 + 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 = None + 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 = None + 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. From 9504d6decc468bb566cb1ed9cd45e7db7952cf92 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Fri, 6 Nov 2026 08:53:54 +0000 Subject: [PATCH 4/8] MCPLGen: fix clang -Winconsistent-missing-override error, flush messages Do not mark the new methods as override, since the other methods of the class are not (clang otherwise gives an error with -Werror). Also flush the output of the messages about the weights, so they are not lost when multiprocessing. Assisted-by: Claude --- .../pkgs/MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) 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 af46f6e..15a30ef 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 @@ -31,9 +31,9 @@ class MCPLGen : public G4Interfaces::ParticleGenBase, void gen(G4Event*); //For MCPLWriter (see G4MCPL/InitialRayCountInfo.hh): - double inputInitialRayCount() const override; - double inputWeightFactor() const override; - bool usedAllInputParticles() const override; + virtual double inputInitialRayCount() const; + virtual double inputWeightFactor() const; + virtual bool usedAllInputParticles() const; virtual bool unlimited() const { return false; }//only has the events available in the input @@ -258,9 +258,11 @@ void MCPLGen::init_weight_scale() //Only print in one process when multiprocessing: const bool doprint = FrameworkGlobals::isParent(); if (!getParameterBoolean("normalise_by_initial_ray_count")) { - if (doprint && nrays != -2.0) + 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) @@ -279,9 +281,11 @@ void MCPLGen::init_weight_scale() throw std::runtime_error("MCPLGen: Can not normalise weights by stat:sum:initial_ray_count"); } m_weight_scale = inputWeightFactor(); - if (doprint) + 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) From 45f8611742e8fd83e2cf4cad843c76d908f9ba6a Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Fri, 6 Nov 2026 13:02:30 +0000 Subject: [PATCH 5/8] ExprParser: define exception destructors out of line The exception classes of ExprParser were header-only, so each module using them had its own copy of their type information. On platforms where libc++ compares type information by address (Intel macOS), exceptions thrown in one module (e.g. for invalid filter expressions) could then not be caught by type in another (e.g. a Python module loaded with RTLD_LOCAL). Defining their destructors in Exception.cc places the type information in the ExprParser library only. Also make testmcplstatsums robust against messages from multiple processes being interleaved in the output. Assisted-by: Claude --- .../pkgs/Utils/ExprParser/libinc/Exception.hh | 16 ++++++++++++++-- .../pkgs/Utils/ExprParser/libsrc/Exception.cc | 8 ++++++++ .../G4MCPLTests/scripts/testmcplstatsums | 11 ++++++++--- 3 files changed, 30 insertions(+), 5 deletions(-) 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 db5c682..d4033ad 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 811ee32..7a61971 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/scripts/testmcplstatsums b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums index e2c52c3..61d83df 100755 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums @@ -28,10 +28,15 @@ def sim( title, infile, *args, expect_ok = True ): cmd = [ 'sb_g4mcpltests_mcplstatsums', f'input_file={infile}', f'--mcpl=grabsrc to {outfile}' ] + list(args) r = subprocess.run( cmd, capture_output = True, text = True, check = False ) - #Only show the lines of MCPLGen and MCPLWriter (others are not relevant here): + #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(): - if line.startswith( ( 'MCPLGen', 'MCPLWriter: Done merging' ) ): - print( ' ', line ) + for m in ( 'MCPLGen: Dividing', 'MCPLGen: Using', 'MCPLGen ERROR', + 'MCPLWriter: Done merging' ): + 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)' ) From cc15501fb41b957125bbcfc0e2260b0623a18315 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Fri, 6 Nov 2026 13:02:30 +0000 Subject: [PATCH 6/8] CI: add Intel macOS (macos-15-intel) to condatest On Intel macOS, libc++ compares type information by address, which makes it sensitive to issues with polymorphic classes or exceptions shared between Python modules, which are not visible on Linux or arm64 macOS. Assisted-by: Claude --- .github/workflows/condatest.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.github/workflows/condatest.yml b/.github/workflows/condatest.yml index 373b6d2..714012a 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: From a27c51452bb9e763551c264ddbec63e0bf48790e Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Fri, 6 Nov 2026 14:14:28 +0000 Subject: [PATCH 7/8] G4MCPL: global registry of stat:sum values for MCPLWriter output files Replace the G4MCPL/InitialRayCountInfo.hh interface (which MCPLWriter found via dynamic_cast on the generator, which failed on Intel macOS since the header-only class had separate type information in each Python module) with a global registry in G4MCPL/OutputStatSums.hh. Components or user code can add stat:sum entries for all files written by MCPLWriter (with a mode which determines how the values add up when the files of multiple processes are merged), and set a factor for the weights written. MCPLGen uses it to transfer stat:sum:initial_ray_count of its input file. The registration happens when its parameters are validated, before the simulation starts and before any fork when multiprocessing. Python code can add constant values with G4MCPLPlugins.MCPLWriter.addStatSum(key,value). Assisted-by: Claude --- doc/source/mcpl.rst | 14 ++++ .../MCPL/G4MCPL/libinc/InitialRayCountInfo.hh | 36 ---------- .../pkgs/MCPL/G4MCPL/libinc/OutputStatSums.hh | 55 ++++++++++++++ .../pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc | 54 ++++++++++++++ .../MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc | 70 ++++++++++++++---- .../G4MCPLPlugins/pymod__MCPLWriter/mod.cc | 72 ++++++++++--------- .../MCPL/G4MCPLPlugins/python/MCPLWriter.py | 11 ++- .../G4MCPLTests/scripts/mcplstatsums | 6 ++ .../G4MCPLTests/scripts/testmcplstatsums | 14 +++- .../G4MCPLTests/scripts/testmcplstatsums.log | 17 +++++ 10 files changed, 263 insertions(+), 86 deletions(-) delete mode 100644 src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/InitialRayCountInfo.hh create mode 100644 src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/OutputStatSums.hh create mode 100644 src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc diff --git a/doc/source/mcpl.rst b/doc/source/mcpl.rst index e0ceb29..b231d09 100644 --- a/doc/source/mcpl.rst +++ b/doc/source/mcpl.rst @@ -123,6 +123,20 @@ weights as in the original simulation. Since the value is only valid when all particles of the input file were used, it is marked as unavailable when this was not the case (e.g. with ``skip_events`` or a limited number of events). +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 transfers ``stat:sum:initial_ray_count``. + Modify MCPL files ----------------- diff --git a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/InitialRayCountInfo.hh b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/InitialRayCountInfo.hh deleted file mode 100644 index 17ce857..0000000 --- a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/InitialRayCountInfo.hh +++ /dev/null @@ -1,36 +0,0 @@ -#ifndef G4MCPL_InitialRayCountInfo_hh -#define G4MCPL_InitialRayCountInfo_hh - -namespace G4MCPL { - - //Interface implemented by particle generators reading particles from MCPL - //files (i.e. G4MCPLPlugins/MCPLGen), with information about the - //stat:sum:initial_ray_count value of the input file. This is used by - //G4MCPLPlugins/MCPLWriter to transfer the value to its output files, with - //the weights written in the same convention as in the input file (i.e. the - //weights of the input file divided by the value are the weights per initial - //ray). - - class InitialRayCountInfo { - public: - virtual ~InitialRayCountInfo() = default; - - //The stat:sum:initial_ray_count value of the input file: -2 if not - //present, -1 if not available, and otherwise the (positive) value. This - //only reads the header of the input file, so it can be called at any time - //after the generator is configured: - virtual double inputInitialRayCount() const = 0; - - //The factor applied to the weights of the input file when generating - //particles (1/initial_ray_count when normalising the weights, otherwise - //1): - virtual double inputWeightFactor() const = 0; - - //Whether all particles of the input file (in this process, when - //multiprocessing) were used: - virtual bool usedAllInputParticles() const = 0; - }; - -} - -#endif 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 0000000..83ec09b --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/OutputStatSums.hh @@ -0,0 +1,55 @@ +#ifndef G4MCPL_OutputStatSums_hh +#define G4MCPL_OutputStatSums_hh + +#include +#include +#include + +//Global (per process) registry of stat:sum values to add to all MCPL files +//written by G4MCPLPlugins/MCPLWriter, and of a factor to apply to the weights +//written in those files. This is for instance used by G4MCPLPlugins/MCPLGen to +//transfer stat:sum:initial_ray_count of its input file, but it can also be +//used by other components or user code to add custom values. +// +//Since MCPL requires the same header entries in all files to merge them (e.g. +//the files of the processes of multi-process jobs), entries must be added +//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 only the parent process + //writes it (the others write 0). If the value is unavailable + //(-1) in any process, it is written as -1 in all processes + //(which makes the merged value unavailable). + }; + + //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(); + + //Factor applied to the Geant4 weights when writing them (default 1). This is + //for instance used by MCPLGen, so the weights are written in the same + //convention as in its input file, when transferring initial_ray_count: + void setOutputWeightFactor( double ); + double outputWeightFactor(); + +} + +#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 0000000..cd2a0ad --- /dev/null +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc @@ -0,0 +1,54 @@ +#include "G4MCPL/OutputStatSums.hh" +#include +#include + +namespace G4MCPL { + namespace { + std::vector& registry() + { + static std::vector s_registry; + return s_registry; + } + double& weightFactor() + { + static double s_factor = 1.0; + return s_factor; + } + } +} + +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" ) + 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::setOutputWeightFactor( double f ) +{ + if ( !( f > 0.0 ) || !std::isfinite(f) ) + throw std::runtime_error("G4MCPL: Output weight factor must be positive and finite"); + weightFactor() = f; +} + +double G4MCPL::outputWeightFactor() +{ + return weightFactor(); +} 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 15a30ef..444c2c4 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,5 +1,5 @@ #include "G4MCPL/G4MCPLUserFlags.hh" -#include "G4MCPL/InitialRayCountInfo.hh" +#include "G4MCPL/OutputStatSums.hh" #include "G4Interfaces/ParticleGenToPy.hh" #include "G4Interfaces/ParticleGenBase.hh" #include "G4Interfaces/FrameworkGlobals.hh" @@ -16,25 +16,20 @@ #include "CLHEP/Geometry/Vector3D.h" #include "CLHEP/Geometry/Transform3D.h" #include +#include #include #include #include #include "G4ParticleTable.hh" #include "G4IonTable.hh" -class MCPLGen : public G4Interfaces::ParticleGenBase, - public G4MCPL::InitialRayCountInfo { +class MCPLGen : public G4Interfaces::ParticleGenBase { public: MCPLGen(); virtual ~MCPLGen(); void init(); void gen(G4Event*); - //For MCPLWriter (see G4MCPL/InitialRayCountInfo.hh): - virtual double inputInitialRayCount() const; - virtual double inputWeightFactor() const; - virtual bool usedAllInputParticles() const; - virtual bool unlimited() const { return false; }//only has the events available in the input protected: @@ -45,7 +40,15 @@ class MCPLGen : public G4Interfaces::ParticleGenBase, 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; + void signalEndOfInput(bool including_current_event); mcpl_file_t m_mcplfile; const mcpl_particle_t * m_p;//particle to be used for next generation bool m_unfiltered; @@ -63,6 +66,11 @@ class MCPLGen : public G4Interfaces::ParticleGenBase, 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; + //Whether all particles of the input file were used (in this process). Shared + //with the function registered in register_output_stat_sums(), which might be + //evaluated after this generator is deleted: + std::shared_ptr m_used_all_input_particles; }; MCPLGen::MCPLGen() @@ -75,7 +83,9 @@ MCPLGen::MCPLGen() m_has_polarisation(false), m_allow_zero_weight(false), m_weight_scale(1.0), - m_input_irc(std::nan("")) + m_input_irc(std::nan("")), + m_registered_output_stat_sums(false), + m_used_all_input_particles(std::make_shared(false)) { m_mcplfile.internal = 0; @@ -144,9 +154,40 @@ 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() +{ + //Transfer stat:sum:initial_ray_count of the input file to MCPL files written + //by MCPLWriter, with the weights in the same convention as in the input file + //(i.e. undoing any normalisation by this generator), so reading the output + //files again gives the same weights. This happens when the parameters are + //validated, which is before the simulation starts (and before any fork when + //multiprocessing), as required by G4MCPL/OutputStatSums.hh: + if (m_registered_output_stat_sums) + return; + m_registered_output_stat_sums = true; + const double nrays = inputInitialRayCount(); + if (nrays == -2.0) + return;//not present + //The value is only valid if all particles of the input file were used: + auto used_all = m_used_all_input_particles; + G4MCPL::addOutputStatSum( "initial_ray_count", + [used_all,nrays]() { return ( nrays > 0.0 && *used_all ) ? nrays : -1.0; }, + G4MCPL::StatSumMode::WholeJob ); + G4MCPL::setOutputWeightFactor( 1.0 / inputWeightFactor() ); +} + void MCPLGen::init() { //In case of multiprocessing we haven't forked yet and each process needs it's @@ -217,9 +258,12 @@ double MCPLGen::inputWeightFactor() const return 1.0; } -bool MCPLGen::usedAllInputParticles() const +void MCPLGen::signalEndOfInput(bool including_current_event) { - return getParameterInt("skip_events") == 0 && reachedLimit(); + //All particles of the input file were used (in this process), unless some + //were skipped: + *m_used_all_input_particles = ( getParameterInt("skip_events") == 0 ); + signalEndOfEvents(including_current_event); } void MCPLGen::delayed_init() @@ -332,7 +376,7 @@ void MCPLGen::gen(G4Event* evt) //Oups, we don't have a particle to generate in this event! Should only //happen if number of (filtered) particles in the file is less than the //number of processes: - signalEndOfEvents(true); + signalEndOfInput(true); return; } @@ -406,7 +450,7 @@ void MCPLGen::gen(G4Event* evt) } } if (!m_p) - signalEndOfEvents(false);//signal that this will be the last event + signalEndOfInput(false);//signal that this will be the last event } 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 3bc1221..d793ce2 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,7 +2,7 @@ #include "G4Interfaces/ComponentHandles.hh" #include "G4MCPL/G4MCPLUserFlags.hh" -#include "G4MCPL/InitialRayCountInfo.hh" +#include "G4MCPL/OutputStatSums.hh" #include "G4Interfaces/FrameworkGlobals.hh" #include "G4Interfaces/GeoConstructBase.hh" #include "G4Interfaces/ParticleGenBase.hh" @@ -21,6 +21,7 @@ #include "G4RegionStore.hh" #include +#include #include #include #include @@ -124,8 +125,7 @@ namespace G4MCPLWriter { EP::Evaluator filter_evaluator, EP::Evaluator flag_evaluator, std::vector> comments_and_blobs, - EventStatsPtr evtstats, - const G4MCPL::InitialRayCountInfo* ircinfo) + EventStatsPtr evtstats) : G4VSensitiveDetector(name), m_expr_builder(builder), m_eval_filter(filter_evaluator), @@ -143,8 +143,6 @@ namespace G4MCPLWriter { m_om(0), m_comments_and_blobs(comments_and_blobs), m_evtstats(evtstats), - m_ircinfo(ircinfo), - m_transfer_irc(false), m_weight_factor(1.0) { std::memset(&m_p,0,sizeof(m_p)); @@ -186,22 +184,25 @@ namespace G4MCPLWriter { mcpl_hdr_add_stat_sum(m_f,"geant4_nevts",(double)m_evtstats->nevts()); mcpl_hdr_add_stat_sum(m_f,"initial_geant4_weights",m_evtstats->sumWeights()); } - if (m_transfer_irc) { - //The initial_ray_count of the input file is only valid if all of its - //particles were used. When multiprocessing, the value is only set in - //the file of the parent process (since the values are added up when - //the files are merged), and the files of all processes must have used - //all their particles (-1 in any of the files makes the merged value - //unavailable): - // - //NB: The value would preferably be set only after merging the files, - //but MCPL currently does not support editing stat:sum values of - //existing files. - const double irc = m_ircinfo->inputInitialRayCount(); - double value = -1.0; - if ( irc > 0.0 && m_ircinfo->usedAllInputParticles() ) - value = FrameworkGlobals::isChild() ? 0.0 : irc; - mcpl_hdr_add_stat_sum(m_f,"initial_ray_count",value); + for ( auto& e : m_statsums ) { + //Final values of the entries from G4MCPL/OutputStatSums.hh (reserved + //in initmcpl()): + double value = e.value(); + if ( !( value >= 0.0 ) || !std::isfinite(value) ) + value = -1.0;//not available + if ( e.mode == G4MCPL::StatSumMode::WholeJob && FrameworkGlobals::isChild() && value != -1.0 ) { + //When multiprocessing, values describing the whole job are only set + //in the file of the parent process, since the values are added up + //when the files are merged (an unavailable value in any of the files + //makes the merged value unavailable). + // + //NB: The value would preferably be set only after merging the files, + //but MCPL does not currently support updating stat:sum values of + //existing files (see mctools/mcpl#254 and + //mctools/simplebuild-dgcode#63). + value = 0.0; + } + mcpl_hdr_add_stat_sum(m_f,e.key.c_str(),value); } if (!FrameworkGlobals::isForked()) mcpl_closeandgzip_outfile(m_f); @@ -353,15 +354,15 @@ namespace G4MCPLWriter { mcpl_hdr_add_stat_sum(m_f,"initial_geant4_weights",-1.0); } - //Transfer the stat:sum:initial_ray_count of an input file (read by - //MCPLGen), with the weights in the same convention as in the input file - //(i.e. undoing any normalisation by MCPLGen), so the weights per initial - //ray are the same when reading the file again: - if ( m_ircinfo && m_ircinfo->inputInitialRayCount() != -2.0 ) { - m_transfer_irc = true; - m_weight_factor = 1.0 / m_ircinfo->inputWeightFactor(); - mcpl_hdr_add_stat_sum(m_f,"initial_ray_count",-1.0); - } + //Values added by other components or user code (e.g. MCPLGen adds the + //stat:sum:initial_ray_count of its input file, with a factor for the + //weights, so they are written in the same convention as in the input + //file). The space for the values is reserved here, and the values are + //set when the file is closed: + m_statsums = G4MCPL::outputStatSums(); + for ( auto& e : m_statsums ) + mcpl_hdr_add_stat_sum(m_f,e.key.c_str(),-1.0); + m_weight_factor = G4MCPL::outputWeightFactor(); if (m_opt_writedoubleprec) mcpl_enable_doubleprec(m_f); if (m_opt_writepolarisation) mcpl_enable_polarisation(m_f); @@ -386,8 +387,7 @@ namespace G4MCPLWriter { MCPLOutputMerger * m_om; std::vector> m_comments_and_blobs; EventStatsPtr m_evtstats; - const G4MCPL::InitialRayCountInfo* m_ircinfo; - bool m_transfer_irc; + std::vector m_statsums; double m_weight_factor; }; @@ -696,8 +696,7 @@ namespace G4MCPLWriter { m_eval_filter, m_eval_flags, createCommentsAndBlobs(), - m_evtstats, - dynamic_cast(m_gen)); + m_evtstats); G4SDManager::GetSDMpointer()->AddNewDetector(m_sd); } auto currentsd = lv->GetSensitiveDetector(); @@ -759,4 +758,9 @@ SBLDPY_MODULE( _MCPLWriter ) .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 ead7478..c514525 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. diff --git a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/mcplstatsums b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/mcplstatsums index 9bc9e74..4a28af9 100755 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/mcplstatsums +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/mcplstatsums @@ -11,6 +11,12 @@ 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") diff --git a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums index 61d83df..0933ecd 100755 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums @@ -11,6 +11,7 @@ # * MCPLWriter transfers stat:sum:initial_ray_count of the input file to the # output file (marked as unavailable unless all particles of the input file # were used), with the weights in the same convention as in the input file. +# * 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. @@ -21,13 +22,16 @@ 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 ): +def sim( title, infile, *args, expect_ok = True, env_extra = None ): print() print( f'==> {title}' ) outfile = 'out.mcpl' cmd = [ 'sb_g4mcpltests_mcplstatsums', f'input_file={infile}', f'--mcpl=grabsrc to {outfile}' ] + list(args) - r = subprocess.run( cmd, capture_output = True, text = True, check = False ) + 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: @@ -76,6 +80,12 @@ sim( 'File with initial_ray_count, for reading again', 'rays.mcpl' ) subprocess.run( [ 'mv', 'out.mcpl.gz', 'again.mcpl.gz' ], check = True ) sim( 'Reading the output again', 'again.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: 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 index 546ea8a..78f5441 100644 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log @@ -77,6 +77,23 @@ Output stat:sum:initial_geant4_weights = 5 Output stat:sum:initial_ray_count = 1000 +==> 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: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5 + Output stat:sum:initial_ray_count = 1000 + Output stat:sum:my_beam_time = 17.5 + +==> 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: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + Output stat:sum:geant4_nevts = 5 + Output stat:sum:initial_geant4_weights = 5 + Output stat:sum:initial_ray_count = 1000 + Output stat:sum:my_beam_time = 17.5 + ==> Filtering a file Output particles: 3 Output stat:sum:initial_ray_count = 1000 From 6af8635d04272105d0400c1e3c069f00e5b54418 Mon Sep 17 00:00:00 2001 From: Thomas Kittelmann Date: Sat, 7 Nov 2026 15:42:23 +0000 Subject: [PATCH 8/8] MCPLWriter: do not claim the initial_ray_count convention by default The stat:sum:initial_ray_count convention of instrument simulations assumes that each particle in a file is an independent ray, which is usually not the case for files written by Geant4. MCPLWriter therefore no longer transfers the initial_ray_count of an MCPLGen input file automatically: - By default, the output files have the weights as used in Geant4, and the value of the input file is recorded as stat:sum:original_initial_ray_count. - With the new opt_irc keyword of --mcpl (or setWriteInitialRayCount(True)), the value is recorded as stat:sum:initial_ray_count, with the weights in the same convention as in the input file (an error if the input file has no valid value). - stat:sum:initial_geant4_weights is in the same convention as the weights in the file. - Values describing the whole job (the above, and custom values added with addStatSum) are divided between the processes when multiprocessing, so the merged file gets the value, rather than being written by the parent process only. The check of whether all particles of the input file were used is removed. Assisted-by: Claude --- doc/source/mcpl.rst | 38 +++++--- .../pkgs/G4/G4Launcher/python/_launcher.py | 8 +- .../pkgs/MCPL/G4MCPL/libinc/OutputStatSums.hh | 39 ++++---- .../pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc | 36 +++++--- .../MCPL/G4MCPLPlugins/pymod__MCPLGen/mod.cc | 37 ++------ .../G4MCPLPlugins/pymod__MCPLWriter/mod.cc | 89 +++++++++++++------ .../MCPL/G4MCPLPlugins/python/MCPLWriter.py | 2 +- .../G4MCPLTests/scripts/testmcplstatsums | 38 +++++--- .../G4MCPLTests/scripts/testmcplstatsums.log | 71 ++++++++++----- .../G4MCPLTests/scripts/testmcplwrite.log | 8 +- 10 files changed, 232 insertions(+), 134 deletions(-) diff --git a/doc/source/mcpl.rst b/doc/source/mcpl.rst index b231d09..ac00af4 100644 --- a/doc/source/mcpl.rst +++ b/doc/source/mcpl.rst @@ -106,22 +106,35 @@ options, by using the module. See an example :sbpkg:`here ` for how to use it. -The MCPL files written by Geant4 contain two ``stat:sum`` values, which describe -the simulation as a whole (independently of how many particles ended up in the +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. - -If the input file read with MCPLGen has a ``stat:sum:initial_ray_count`` value, -it is transferred to the output file, and the weights are written in the same -convention as in the input file (i.e. multiplied by the value if MCPLGen divided -the weights by it). So reading the output file with MCPLGen again gives the same -weights as in the original simulation. Since the value is only valid when all -particles of the input file were used, it is marked as unavailable when this -was not the case (e.g. with ``skip_events`` or a limited number of events). + 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 @@ -135,7 +148,8 @@ the :ref:`sim script `, before the simulation is started: 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 transfers ``stat:sum:initial_ray_count``. +which is also how MCPLGen provides the ``stat:sum:initial_ray_count`` of its +input file. Modify MCPL files ----------------- 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 95a7110..6f57a55 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 index 83ec09b..fb2f637 100644 --- a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/OutputStatSums.hh +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libinc/OutputStatSums.hh @@ -5,17 +5,18 @@ #include #include -//Global (per process) registry of stat:sum values to add to all MCPL files -//written by G4MCPLPlugins/MCPLWriter, and of a factor to apply to the weights -//written in those files. This is for instance used by G4MCPLPlugins/MCPLGen to -//transfer stat:sum:initial_ray_count of its input file, but it can also be -//used by other components or user code to add custom values. +//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), entries must be added -//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). +//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 { @@ -23,10 +24,9 @@ namespace G4MCPL { //(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 only the parent process - //writes it (the others write 0). If the value is unavailable - //(-1) in any process, it is written as -1 in all processes - //(which makes the merged value unavailable). + 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 @@ -44,11 +44,14 @@ namespace G4MCPL { }; const std::vector& outputStatSums(); - //Factor applied to the Geant4 weights when writing them (default 1). This is - //for instance used by MCPLGen, so the weights are written in the same - //convention as in its input file, when transferring initial_ray_count: - void setOutputWeightFactor( double ); - double outputWeightFactor(); + //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(); } diff --git a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc index cd2a0ad..e5368a0 100644 --- a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPL/libsrc/OutputStatSums.cc @@ -9,10 +9,15 @@ namespace G4MCPL { static std::vector s_registry; return s_registry; } - double& weightFactor() + struct InputInfo { + bool present = false; + double initial_ray_count = -1.0; + double weight_factor = 1.0; + }; + InputInfo& inputInfo() { - static double s_factor = 1.0; - return s_factor; + static InputInfo s_info; + return s_info; } } } @@ -26,7 +31,8 @@ void G4MCPL::addOutputStatSum( const std::string& key, 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" ) + 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 ) @@ -41,14 +47,20 @@ const std::vector& G4MCPL::outputStatSums() return registry(); } -void G4MCPL::setOutputWeightFactor( double f ) +void G4MCPL::setInputInitialRayCount( double initial_ray_count, double weight_factor ) { - if ( !( f > 0.0 ) || !std::isfinite(f) ) - throw std::runtime_error("G4MCPL: Output weight factor must be positive and finite"); - weightFactor() = f; + 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; } -double G4MCPL::outputWeightFactor() -{ - return weightFactor(); -} +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 444c2c4..671946b 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 @@ -16,7 +16,6 @@ #include "CLHEP/Geometry/Vector3D.h" #include "CLHEP/Geometry/Transform3D.h" #include -#include #include #include #include @@ -48,7 +47,7 @@ class MCPLGen : public G4Interfaces::ParticleGenBase { //(in this process, when multiprocessing): double inputInitialRayCount() const; double inputWeightFactor() const; - void signalEndOfInput(bool including_current_event); + mcpl_file_t m_mcplfile; const mcpl_particle_t * m_p;//particle to be used for next generation bool m_unfiltered; @@ -67,10 +66,6 @@ class MCPLGen : public G4Interfaces::ParticleGenBase { 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; - //Whether all particles of the input file were used (in this process). Shared - //with the function registered in register_output_stat_sums(), which might be - //evaluated after this generator is deleted: - std::shared_ptr m_used_all_input_particles; }; MCPLGen::MCPLGen() @@ -84,8 +79,7 @@ MCPLGen::MCPLGen() m_allow_zero_weight(false), m_weight_scale(1.0), m_input_irc(std::nan("")), - m_registered_output_stat_sums(false), - m_used_all_input_particles(std::make_shared(false)) + m_registered_output_stat_sums(false) { m_mcplfile.internal = 0; @@ -168,24 +162,18 @@ bool MCPLGen::validateParameters() void MCPLGen::register_output_stat_sums() { - //Transfer stat:sum:initial_ray_count of the input file to MCPL files written - //by MCPLWriter, with the weights in the same convention as in the input file - //(i.e. undoing any normalisation by this generator), so reading the output - //files again gives the same weights. This happens when the parameters are + //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), as required by G4MCPL/OutputStatSums.hh: + //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 - //The value is only valid if all particles of the input file were used: - auto used_all = m_used_all_input_particles; - G4MCPL::addOutputStatSum( "initial_ray_count", - [used_all,nrays]() { return ( nrays > 0.0 && *used_all ) ? nrays : -1.0; }, - G4MCPL::StatSumMode::WholeJob ); - G4MCPL::setOutputWeightFactor( 1.0 / inputWeightFactor() ); + G4MCPL::setInputInitialRayCount( nrays, inputWeightFactor() ); } void MCPLGen::init() @@ -258,13 +246,6 @@ double MCPLGen::inputWeightFactor() const return 1.0; } -void MCPLGen::signalEndOfInput(bool including_current_event) -{ - //All particles of the input file were used (in this process), unless some - //were skipped: - *m_used_all_input_particles = ( getParameterInt("skip_events") == 0 ); - signalEndOfEvents(including_current_event); -} void MCPLGen::delayed_init() { @@ -376,7 +357,7 @@ void MCPLGen::gen(G4Event* evt) //Oups, we don't have a particle to generate in this event! Should only //happen if number of (filtered) particles in the file is less than the //number of processes: - signalEndOfInput(true); + signalEndOfEvents(true); return; } @@ -450,7 +431,7 @@ void MCPLGen::gen(G4Event* evt) } } if (!m_p) - signalEndOfInput(false);//signal that this will be the last event + signalEndOfEvents(false);//signal that this will be the last event } 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 d793ce2..af7dc4f 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 @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -119,6 +120,7 @@ namespace G4MCPLWriter { bool writepolarisation, bool writeuserflags, double universalweight, + bool writeinitialraycount, bool writeonvolexit, KillStrategy killstrategy, BuilderPtr builder, @@ -143,6 +145,7 @@ namespace G4MCPLWriter { m_om(0), 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)); @@ -179,30 +182,20 @@ namespace G4MCPLWriter { return; m_closed = true; assert(m_f.internal); + //Final values of the stat:sum entries reserved in initmcpl(): if (m_evtstats) { - //Final values of the stat:sum entries reserved in initmcpl(): mcpl_hdr_add_stat_sum(m_f,"geant4_nevts",(double)m_evtstats->nevts()); - mcpl_hdr_add_stat_sum(m_f,"initial_geant4_weights",m_evtstats->sumWeights()); + //(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 ) { - //Final values of the entries from G4MCPL/OutputStatSums.hh (reserved - //in initmcpl()): - double value = e.value(); - if ( !( value >= 0.0 ) || !std::isfinite(value) ) - value = -1.0;//not available - if ( e.mode == G4MCPL::StatSumMode::WholeJob && FrameworkGlobals::isChild() && value != -1.0 ) { - //When multiprocessing, values describing the whole job are only set - //in the file of the parent process, since the values are added up - //when the files are merged (an unavailable value in any of the files - //makes the merged value unavailable). - // - //NB: The value would preferably be set only after merging the files, - //but MCPL does not currently support updating stat:sum values of - //existing files (see mctools/mcpl#254 and - //mctools/simplebuild-dgcode#63). - value = 0.0; - } - mcpl_hdr_add_stat_sum(m_f,e.key.c_str(),value); + 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); @@ -354,15 +347,28 @@ namespace G4MCPLWriter { mcpl_hdr_add_stat_sum(m_f,"initial_geant4_weights",-1.0); } - //Values added by other components or user code (e.g. MCPLGen adds the - //stat:sum:initial_ray_count of its input file, with a factor for the - //weights, so they are written in the same convention as in the input - //file). The space for the values is reserved here, and the values are - //set when the file is closed: + //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); - m_weight_factor = G4MCPL::outputWeightFactor(); if (m_opt_writedoubleprec) mcpl_enable_doubleprec(m_f); if (m_opt_writepolarisation) mcpl_enable_polarisation(m_f); @@ -388,7 +394,23 @@ namespace G4MCPLWriter { 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() { @@ -420,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; @@ -431,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; } @@ -444,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) @@ -477,6 +502,12 @@ namespace G4MCPLWriter { 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) @@ -642,7 +673,9 @@ namespace G4MCPLWriter { } 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"); + " 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); @@ -690,6 +723,7 @@ namespace G4MCPLWriter { m_opt_writepolarisation, m_opt_writeuserflags, m_opt_universalweight, + m_opt_writeinitialraycount, m_opt_writeonvolexit, m_opt_killstrategy, m_expr_builder, @@ -755,6 +789,7 @@ 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(); } ); 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 c514525..c950e53 100644 --- a/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/python/MCPLWriter.py +++ b/src/simplebuild_dgcode/data/pkgs/MCPL/G4MCPLPlugins/python/MCPLWriter.py @@ -19,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_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums index 0933ecd..d6bc2ec 100755 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums @@ -8,9 +8,12 @@ # (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 transfers stat:sum:initial_ray_count of the input file to the -# output file (marked as unavailable unless all particles of the input file -# were used), with the weights in the same convention as in the input file. +# * 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. @@ -22,12 +25,12 @@ 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 ): +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=grabsrc to {outfile}' ] + list(args) + f'--mcpl={mcplkw} to {outfile}' ] + list(args) import os env = dict( os.environ ) env.update( env_extra or {} ) @@ -37,7 +40,7 @@ def sim( title, infile, *args, expect_ok = True, env_extra = None ): #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: Done merging', 'MCPLWriter: Writing stat:sum' ): if m in line: print( ' ', line[line.index(m):] ) break @@ -72,13 +75,26 @@ 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 ) -#The initial_ray_count of the input file is transferred to the output file if -#all particles of the input file were used, with the weights in the same -#convention as in the input file. So reading the output again gives the same -#weights in Geant4: -sim( 'File with initial_ray_count, for reading again', 'rays.mcpl' ) +#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): 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 index 78f5441..f6e001e 100644 --- a/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log +++ b/src/simplebuild_dgcode/data/pkgs_val/UnitTests/G4MCPLTests/scripts/testmcplstatsums.log @@ -1,47 +1,47 @@ ==> 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: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + 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:initial_ray_count = 1000 + 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:initial_ray_count = 1000 + 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: [1000.0, 1000.0, 1000.0] + Output weights: [1.0, 1.0, 1.0] Output stat:sum:geant4_nevts = 3 Output stat:sum:initial_geant4_weights = 3 - Output stat:sum:initial_ray_count = + 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: [1000.0, 1000.0] + Output weights: [1.0, 1.0] Output stat:sum:geant4_nevts = 2 Output stat:sum:initial_geant4_weights = 2 - Output stat:sum:initial_ray_count = + 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: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + 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:initial_ray_count = 1000 + 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: [1000.0, 1000.0, 1000.0] + Output weights: [1.0, 1.0, 1.0] Output stat:sum:geant4_nevts = 3 Output stat:sum:initial_geant4_weights = 3 - Output stat:sum:initial_ray_count = + 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] @@ -57,42 +57,73 @@ 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 = + 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, for reading again +==> 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 = 5 + Output stat:sum:initial_geant4_weights = 5000 Output stat:sum:initial_ray_count = 1000 -==> Reading the output again +==> 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 = 5 + Output stat:sum:initial_geant4_weights = 5000 Output stat:sum:initial_ray_count = 1000 -==> Custom stat:sum value +==> 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 = 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: [1000.0, 1000.0, 1000.0, 1000.0, 1000.0] + 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:initial_ray_count = 1000 Output stat:sum:my_beam_time = 17.5 + Output stat:sum:original_initial_ray_count = 1000 ==> Filtering a file Output particles: 3 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 774cf30..51f8ebc 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,7 +871,7 @@ Opened MCPL file target.mcpl.gz: Basic info Format : MCPL-3 No. of particles : 0 - Header storage : 1197 bytes + Header storage : 1325 bytes Data storage : 0 bytes Custom meta data @@ -885,7 +885,7 @@ 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" + -> 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 @@ -906,7 +906,7 @@ Opened MCPL file recordfwd.mcpl.gz: Basic info Format : MCPL-3 No. of particles : 6 - Header storage : 1124 bytes + Header storage : 1252 bytes Data storage : 192 bytes Custom meta data @@ -920,7 +920,7 @@ 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" + -> 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