diff --git a/docs/manual/geombench.md b/docs/manual/geombench.md new file mode 100644 index 00000000..28bddc23 --- /dev/null +++ b/docs/manual/geombench.md @@ -0,0 +1,233 @@ +(manual-geombench)= + +# Geometry benchmarking + +The `remage-geombench` command-line tool helps identify performance bottlenecks +in detector geometries during the design phase. It systematically samples points +on three grids of starting positions and simulates geantinos through the volume. +It returns the median simulation time per geantino. Based on these overview +plots and statistics are generated. + +:::{tip} + +Use this tool during geometry development to identify components that may slow +down simulations. Complex nested structures, boolean operations, and highly +detailed volumes can significantly increase navigation time. + +::: + +## Overview + +The geometry benchmark works by: + +1. Loading a GDML geometry file +2. Creating three 2D grids in XY, XZ, and YZ direction according to grid + specifications +3. Sampling geantinos and estimating the median simulation time per grid point +4. Generating statistics and visualizations showing where navigation is slow + +This information helps you optimize geometries before running full physics +simulations, potentially saving significant computation time. + +## Basic usage + +The simplest usage requires only a GDML file: + +```console +$ remage-geombench detector.gdml +``` + +This will: + +- Benchmark the entire geometry +- Use default settings (10M events, 1 mm grid spacing, 25% buffer) +- Save results to `detector.lh5` in the current directory +- Print summary statistics to the console and generate overview plots in the + current directory + +## Command-line options + +### Geometry selection + +**`geometry`** (required) : Path to the GDML geometry file to benchmark. + +**`--logical-volume NAME`** (optional) : Extract and benchmark only a specific +logical volume including daughters from the geometry. Useful for isolating +performance issues in complex assemblies. + +Example: + +```console +$ remage-geombench l1000.gdml --logical-volume V0101 +``` + +### Grid configuration + +**`--grid-increment SPACING`** (default: `1`) : Uniform spacing between grid +points in millimeters for all dimensions. + +Example (coarse grid for quick tests): + +```console +$ remage-geombench l1000.gdml --grid-increment 5 +``` + +**`--grid-increments DICT`** (optional) : Specify different spacing per +dimension using a Python dictionary literal. Overrides `--grid-increment` if +provided. + +Example: + +```console +$ remage-geombench l1000.gdml --grid-increments "{'x': 1.0, 'y': 2.0, 'z': 0.5}" +``` + +### Buffer region + +**`--buffer-fraction FRACTION`** (default: `0.25`) : Fractional buffer space +around the geometry. A value of 0.25 adds 12.5% extra space on each side, +creating a world volume large enough to contain the geometry with margin. + +Example (tighter bounds): + +```console +$ remage-geombench l1000.gdml --buffer-fraction 0.1 +``` + +### Simulation control + +**`--num-events N`** (default: `10000000`) : Number of navigation events to +simulate. Higher values give more stable statistics but take longer. + +Example (quick test): + +```console +$ remage-geombench l1000.gdml --num-events 1000000 +``` + +**`--output-dir PATH`** (default: `./`) : Directory to store output files. + +Example: + +```console +$ remage-geombench l1000.gdml --output-dir ./benchmark_results/ +``` + +**`--dry-run`** : Generate and display the macro file without running the +simulation. Useful for verifying configuration before long runs. + +Example: + +```console +$ remage-geombench l1000.gdml --dry-run +``` + +## Interpreting results + +The tool generates an LH5 output file containing spatial navigation performance +data and prints summary statistics: + +``` +Geometry Benchmark Analysis Results: +mean_navigation_time: 2.45 ns +median_navigation_time: 1.89 ns +std_navigation_time: 1.23 ns +max_navigation_time: 45.67 ns +min_navigation_time: 0.34 ns +``` + +In addition, the visualization of the data can help identify hotspots. + +### Key metrics + +**Mean navigation time** : Average time spent navigating to each point. Lower is +better. + +**Max navigation time** : Slowest navigation time encountered. High values +indicate geometry hotspots. + +**Standard deviation** : Variability in navigation time. High values suggest +non-uniform complexity. + +:::{warning} + +Navigation times are relative and depend on the hardware. Focus on identifying +spatial patterns and relative differences between geometric components rather +than absolute timing values. + +In addition, these statistics also depend on the ratio of empty space to actual +material. If one is not interested in potential slowdowns far away from the +object, one should choose a small enough buffer to reduce the impact of the +empty space. + +::: + +## Workflow examples + +### Full detector benchmark + +Benchmark an entire detector assembly with default settings: + +```console +$ remage-geombench l1000.gdml --output-dir benchmarks/ +``` + +### Component isolation + +Test a specific problematic component: + +```console +$ remage-geombench l1000.gdml --logical-volume V0101 \ + --grid-increment 0.5 --output-dir component_tests/ +``` + +This generates a `part_{logical-volume}.lh5` and +`part_{logical-volume}_[...].pdf` files in the output directory. + +### Anisotropic sampling + +For elongated geometries, use different grid spacing per dimension: + +```console +$ remage-geombench l1000.gdml \ + --grid-increments "{'x': 0.5, 'y': 0.5, 'z': 2.0}" +``` + +## Performance optimization tips + +Based on benchmark results, consider these optimization strategies: + +1. **Simplify boolean operations**: Multiple nested unions/subtractions are + expensive. Consider alternative representations. + +2. **Reduce tessellated solid complexity**: Decrease facet count where possible + without losing essential features. + +3. **Material boundaries**: Excessive material changes force more boundary + checks. Consolidate materials where physics allows. + +4. **Envelope optimization**: Proper envelope volumes can accelerate navigation + in complex assemblies. + +## Technical details + +The benchmark uses the `GeomBench` generator, which systematically steps through +three 2D grids and measures median simulation time per grid point. The geometry +is automatically wrapped in a world volume sized to contain it with the +specified buffer fraction. + +When extracting a specific logical volume, the tool: + +1. Identifies the volume in the GDML registry +2. Copies all dependent resources (materials, solids, daughter volumes) +3. Creates a minimal world volume containing only the extracted component +4. Applies the buffer and generates a temporary GDML file for benchmarking + +This allows testing individual components without the overhead of the full +geometry. + +## See also + +- {ref}`manual-geometry` – General geometry setup +- {ref}`manual-running` – Running simulations +- {ref}`manual-output` – Output data formats diff --git a/docs/manual/index.md b/docs/manual/index.md index bcfd0397..7950f5b7 100644 --- a/docs/manual/index.md +++ b/docs/manual/index.md @@ -7,6 +7,7 @@ install.md containers.md running.md geometry.md +geombench.md physicslist.md confinement.md generators.md diff --git a/docs/rmg-commands.md b/docs/rmg-commands.md index 7971cb31..92efd396 100644 --- a/docs/rmg-commands.md +++ b/docs/rmg-commands.md @@ -536,7 +536,7 @@ Select event generator * **Parameter** – `generator` * **Parameter type** – `s` * **Omittable** – `False` - * **Candidates** – `G4gun GPS BxDecay0 FromFile CosmicMuons MUSUNCosmicMuons UserDefined Undefined` + * **Candidates** – `G4gun GPS BxDecay0 FromFile CosmicMuons MUSUNCosmicMuons UserDefined GeomBench Undefined` * **Allowed states** – `PreInit Idle` ## `/RMG/Generator/Confinement/` @@ -1266,6 +1266,7 @@ Commands for controlling the simulation output * `/RMG/Output/IsotopeFilter/` – Commands for filtering event out by created isotopes. * `/RMG/Output/Track/` – Commands for controlling output of track vertices. * `/RMG/Output/ParticleFilter/` – Commands for filtering particles out by PDG encoding. +* `/RMG/Output/Benchmark/` – Commands for controlling geometry navigation benchmark output. **Commands:** @@ -2072,6 +2073,10 @@ Add a physics process by name. This will only keep the specified particles when * **Omittable** – `False` * **Allowed states** – `Idle` +## `/RMG/Output/Benchmark/` + +Commands for controlling geometry navigation benchmark output. + ## `/RMG/GrabmayrGammaCascades/` Control Peters gamma cascade model diff --git a/include/RMGGeomBench.hh b/include/RMGGeomBench.hh new file mode 100644 index 00000000..4e3e7a8e --- /dev/null +++ b/include/RMGGeomBench.hh @@ -0,0 +1,85 @@ +// Copyright (C) 2025 Moritz Neuberger +#ifndef _RMG_GEOMBENCH_HH_ +#define _RMG_GEOMBENCH_HH_ + +#include +#include +#include +#include + +#include "CLHEP/Units/SystemOfUnits.h" +#include "G4AnalysisManager.hh" +#include "G4GenericMessenger.hh" +#include "G4LogicalVolume.hh" +#include "G4ParticleGun.hh" +#include "G4VPhysicalVolume.hh" + +#include "RMGGeomBenchOutputScheme.hh" +#include "RMGVGenerator.hh" + +namespace u = CLHEP; + +class RMGGeomBench : public RMGVGenerator { + + public: + + RMGGeomBench(); + ~RMGGeomBench(); + + RMGGeomBench(RMGGeomBench const&) = delete; + RMGGeomBench& operator=(RMGGeomBench const&) = delete; + RMGGeomBench(RMGGeomBench&&) = delete; + RMGGeomBench& operator=(RMGGeomBench&&) = delete; + + void GeneratePrimaries(G4Event* event) override; + void SetParticlePosition(G4ThreeVector) override{}; + void RecordBatchTime(size_t pixel_idx, double batch_time); + void SaveAllPixels(); + + void BeginOfRunAction(const G4Run* r) override; + void EndOfRunAction(const G4Run* r) override; + + private: + + // Helper to find the benchmark output scheme if it's active + RMGGeomBenchOutputScheme* GetBenchmarkOutputScheme(); + + std::unique_ptr fGun = nullptr; + + std::unique_ptr fMessenger = nullptr; + void DefineCommands(); + + long totalnevents; + size_t totalnpixels; + int npixelsperrow; + int neventsperpixel; + double cubesize; + + // Configurable sampling parameters (user-specified increments) + G4ThreeVector user_increment; + G4ThreeVector sampling_width; + + // Calculated number of pixels based on increments and widths + size_t npixels_x; + size_t npixels_y; + size_t npixels_z; + size_t ID; + + double starttime; + double currenttime; + double bunchstarttime; + + // For tracking batches and calculating median + std::vector> pixel_batch_times; // One vector of batch times per pixel + int events_per_bunch; + int total_batch_rounds; + int current_batch_event; + int current_pixel_index; + int current_batch_round; + + G4ThreeVector origin; + G4ThreeVector limit; + G4ThreeVector increment; +}; + +#endif diff --git a/include/RMGGeomBenchOutputScheme.hh b/include/RMGGeomBenchOutputScheme.hh new file mode 100644 index 00000000..d817deb9 --- /dev/null +++ b/include/RMGGeomBenchOutputScheme.hh @@ -0,0 +1,70 @@ +// Copyright (C) 2025 Moritz Neuberger +// +// This program is free software: you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with this program. If not, see . + +#ifndef _RMG_GEOMBENCH_OUTPUT_SCHEME_HH_ +#define _RMG_GEOMBENCH_OUTPUT_SCHEME_HH_ + +#include + +#include "G4AnalysisManager.hh" +#include "G4GenericMessenger.hh" + +#include "RMGVOutputScheme.hh" + +class G4Event; + +/** @brief Output scheme for geometry navigation benchmark data. + * + * @details This output scheme records timing information for geometry navigation + * benchmarking across three orthogonal planes (XZ, YZ, XY). For each pixel in the + * benchmark grid, it stores: + * - Position coordinates (X, Y, Z depending on the plane) + * - Navigation time in seconds + * + * The benchmark data is stored in three separate auxiliary ntuples, one for each plane. + */ +class RMGGeomBenchOutputScheme : public RMGVOutputScheme { + + public: + + RMGGeomBenchOutputScheme(); + + /** @brief Sets the names of the output columns, invoked in @c RMGRunAction::SetupAnalysisManager */ + void AssignOutputNames(G4AnalysisManager* ana_man) override; + + /** @brief Store benchmark pixel data - called when a pixel completes sampling */ + void SavePixel(int plane_id, double x, double y, double z, double time); + + /** @brief Always store benchmark data regardless of event filtering */ + [[nodiscard]] bool StoreAlways() const override { return true; } + + protected: + + [[nodiscard]] std::string GetNtupleName(RMGDetectorMetadata) const override { + throw std::logic_error("benchmark output scheme has no detectors"); + } + + private: + + std::unique_ptr fMessenger; + void DefineCommands(); + + // Ntuple IDs for the three benchmark planes (XZ, YZ, XY) + int fNtupleIDs[3] = {-1, -1, -1}; +}; + +#endif + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/include/RMGMasterGenerator.hh b/include/RMGMasterGenerator.hh index a9d431e0..50d1a5be 100644 --- a/include/RMGMasterGenerator.hh +++ b/include/RMGMasterGenerator.hh @@ -49,6 +49,7 @@ class RMGMasterGenerator : public G4VUserPrimaryGeneratorAction { kCosmicMuons, ///< A simple cosmic muon generator. kMUSUNCosmicMuons, ///< The MUSUN-based cosmic muon generator. kUserDefined, ///< A user-specified custom generator. + kGeomBench, ///< The benchmark generator. kUndefined ///< Undefined generator mode. }; diff --git a/include/RMGUserInit.hh b/include/RMGUserInit.hh index 93a187a6..7d3f982d 100644 --- a/include/RMGUserInit.hh +++ b/include/RMGUserInit.hh @@ -26,6 +26,7 @@ #include "G4UserSteppingAction.hh" #include "G4UserTrackingAction.hh" +#include "RMGGeomBenchOutputScheme.hh" #include "RMGGeometryCheckOutputScheme.hh" #include "RMGIsotopeFilterScheme.hh" #include "RMGLog.hh" @@ -155,6 +156,7 @@ class RMGUserInit { AddOptionalOutputScheme("ParticleFilter"); AddOptionalOutputScheme("Track"); AddOptionalOutputScheme("GeometryCheck"); + AddOptionalOutputScheme("GeomBench"); } /** diff --git a/pyproject.toml b/pyproject.toml index 3945de58..caa7b67c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -38,6 +38,8 @@ dependencies = [ "pygama >=2.2.3", "colorlog", "rich", + "mplhep", + "hist" ] [tool.hatch.metadata] @@ -71,6 +73,7 @@ docs = [ [project.scripts] remage = "remage.cli:remage_cli" +remage-geombench = "remage.geombench.cli:remage_geombench_cli" [project.urls] Homepage = "https://github.com/legend-exp/remage" diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 64495e81..15ea76f1 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -13,7 +13,12 @@ set(PYTHON_SOURCES ${_r}/python/remage/cli.py ${_r}/python/remage/ipc.py ${_r}/python/remage/logging.py - ${_r}/python/remage/post_proc.py) + ${_r}/python/remage/post_proc.py + ${_r}/python/remage/geombench/__init__.py + ${_r}/python/remage/geombench/__main__.py + ${_r}/python/remage/geombench/cli.py + ${_r}/python/remage/geombench/gdml_handling.py + ${_r}/python/remage/geombench/summary_generator.py) # cmake-format: on option(RMG_CONDA_BUILD "Adapt for building on conda-forge." OFF) diff --git a/python/install-remage-python.sh b/python/install-remage-python.sh index 0353617e..05f12f96 100644 --- a/python/install-remage-python.sh +++ b/python/install-remage-python.sh @@ -30,6 +30,7 @@ if [[ "$RMG_PYTHON_PKG_MANAGER" == "uv" ]]; then # install our package into this venv python -m uv -q pip install --reinstall "${CMAKE_SOURCE_DIR}" ln -fs "${INSTALL_VENV_DIR}/bin/remage" "${CMAKE_INSTALL_PREFIX}/bin/remage" + ln -fs "${INSTALL_VENV_DIR}/bin/remage-geombench" "${CMAKE_INSTALL_PREFIX}/bin/remage-geombench" # create and a user-installable wheel python -m uv build --wheel "${CMAKE_SOURCE_DIR}" -o "${CMAKE_INSTALL_PREFIX}/share" diff --git a/python/remage/geombench/__init__.py b/python/remage/geombench/__init__.py new file mode 100644 index 00000000..9cc7dd20 --- /dev/null +++ b/python/remage/geombench/__init__.py @@ -0,0 +1,17 @@ +# Copyright (C) 2025 Moritz Neuberger +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see . + + +"""Geometry benchmark utilities for remage.""" diff --git a/python/remage/geombench/__main__.py b/python/remage/geombench/__main__.py new file mode 100644 index 00000000..b5ea3349 --- /dev/null +++ b/python/remage/geombench/__main__.py @@ -0,0 +1,24 @@ +# Copyright (C) 2025 Moritz Neuberger +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see . + + +from __future__ import annotations + +import sys + +from remage.geombench.cli import remage_geombench_cli + +if __name__ == "__main__": + sys.exit(remage_geombench_cli()) diff --git a/python/remage/geombench/cli.py b/python/remage/geombench/cli.py new file mode 100644 index 00000000..453f51c1 --- /dev/null +++ b/python/remage/geombench/cli.py @@ -0,0 +1,228 @@ +# Copyright (C) 2025 Moritz Neuberger +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see . + + +from __future__ import annotations + +import argparse +import ast +from pathlib import Path + +from .. import logging as rmg_logging +from ..cli import remage_run +from .gdml_handling import ( + generate_tmp_gdml_geometry, + load_gdml_geometry, +) +from .summary_generator import SummaryGenerator + + +def generate_macro(args, output_file_stem: str = "") -> str: + """Generate a macro file for the remage geometry benchmark. + + This function creates a macro file based on the provided command-line + arguments, which will be used to configure the remage simulation for the + geometry benchmark. + + Parameters + ---------- + args + Parsed command-line arguments. + Consisting of: + - geometry: Path to the geometry file. + - output_dir: Directory to store benchmark output files. + - num_events: Number of events to simulate. + - grid_increment: Increment between grid points in mm. + - grid_increments: Optional dict of increments per dimension in mm. + - dry_run: If set, only generate the macro file without running remage. + + Returns + ------- + str + Path to the generated macro file. + """ + + geometry_path = Path(args.geometry) + filename = geometry_path.absolute() + output_dir = Path(args.output_dir) + output_file = output_dir / (output_file_stem + ".lh5") + if issubclass(type(args.grid_increments), str) and args.grid_increments: + grid_increments = ast.literal_eval(args.grid_increments) + increment_x = grid_increments.get("x", args.grid_increment) + increment_y = grid_increments.get("y", args.grid_increment) + increment_z = grid_increments.get("z", args.grid_increment) + else: + increment_x = args.grid_increment + increment_y = args.grid_increment + increment_z = args.grid_increment + n_events = args.num_events + + commands = [ + f"/RMG/Geometry/IncludeGDMLFile {filename}", + "/RMG/Output/ActivateOutputScheme GeomBench", + f"/RMG/Output/FileName {output_file}", + "/run/initialize", + "/RMG/Generator/Select GeomBench", + f"/RMG/Generator/Benchmark/IncrementX {increment_x} mm", + f"/RMG/Generator/Benchmark/IncrementY {increment_y} mm", + f"/RMG/Generator/Benchmark/IncrementZ {increment_z} mm", + f"/run/beamOn {n_events}", + ] + macro_content = "\n".join(commands) + + if args.dry_run: + logger = rmg_logging.setup_log() + logger.info("Dry run enabled. Generated macro content:") + logger.info(macro_content) + return "" + + return macro_content + + +def remage_geombench_cli(external_args: list[str] | None = None) -> int: + """Command-line interface for the remage geometry benchmark. + + This function parses command-line arguments, sets up the macros, runs + the geometry benchmark using the remage simulation framework, and analyzes + the output plus provides a summary of the results. + + Returns + ------- + int + Exit code of the program. Returns 0 on success, non-zero on failure. + """ + parser = argparse.ArgumentParser( + description="Run remage geometry benchmark and analyze results.", + allow_abbrev=False, + ) + parser.add_argument( + "geometry", + type=str, + help="Path to the geometry file to be used in the benchmark.", + ) + parser.add_argument( + "--logical-volume", + type=str, + default="", + help="In case one is interested in a specific logical volume in a complex GDML geometry, only this volume will be benchmarked.", + ) + parser.add_argument( + "--buffer-fraction", + type=float, + default=0.25, + help="Fractional buffer to add around the geometry. For example, 0.25 adds 12.5%% extra space on each side.", + ) + parser.add_argument( + "--output-dir", + type=str, + default="./", + help="Directory to store benchmark output files.", + ) + parser.add_argument( + "--num-events", + type=int, + default=10000000, + help="Number of events to simulate in the benchmark.", + ) + parser.add_argument( + "--grid-increment", + type=float, + default=1, + help="Increment between grid points in the benchmark geometry given in mm. The same for all dimensions.", + ) + parser.add_argument( + "--grid-increments", + type=str, + default="", + help="Increment of individual grid point distances per dimension given in mm. Example: \"{'x': 1., 'y': 2., 'z': 3}\"", + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="If set, only generate the macro file without running remage.", + ) + + # merge external args if provided + args = parser.parse_args(args=external_args) + + logger = rmg_logging.setup_log() + + tmp_gdml_file = "" + + original_gdml_dict = load_gdml_geometry(Path(args.geometry)) + output_file_stem = Path(args.geometry).stem + + # Extract specific component if requested, otherwise use full geometry + if args.logical_volume != "": + registry = original_gdml_dict["registry"] + if args.logical_volume not in registry.logicalVolumeDict: + msg = f"Logical volume '{args.logical_volume}' not found in the geometry registry." + raise ValueError(msg) + + geometry_to_benchmark = { + "object_lv": registry.logicalVolumeDict[args.logical_volume], + "registry": registry, + } + + object_name = f"{args.logical_volume}_extracted" + output_file_stem = f"part_{args.logical_volume}_{output_file_stem}" + else: + geometry_to_benchmark = original_gdml_dict + object_name = "object_lv" + + # Generate temporary GDML with buffered world volume + tmp_gdml_file = generate_tmp_gdml_geometry( + geometry_to_benchmark, + buffer_fraction=args.buffer_fraction, + object_name=object_name, + ) + args.geometry = str(tmp_gdml_file) + + macro_content = generate_macro(args, output_file_stem=output_file_stem) + + if args.dry_run: + return 0 + + try: + # run remage + ec, _ = remage_run(macros=macro_content) + + sim_output_file = Path(args.output_dir) / (output_file_stem + ".lh5") + + if ec != 0 and not sim_output_file.exists(): + logger.error("Remage simulation failed.") + return int(ec) + + sum_gen = SummaryGenerator( + sim_output_file=sim_output_file, + args=args, + output_file_stem=output_file_stem, + ) + analysis_results = sum_gen.perform_analysis() + logger.info("Geometry Benchmark Analysis Results:") + for key, value in analysis_results.items(): + msg = f"{key}: {value}" + logger.info(msg) + + except Exception as e: + msg = f"An error occurred: {e}" + logger.error(msg) + return 1 + + finally: + if tmp_gdml_file: + Path(tmp_gdml_file).unlink(missing_ok=True) + + return 0 diff --git a/python/remage/geombench/gdml_handling.py b/python/remage/geombench/gdml_handling.py new file mode 100644 index 00000000..3a44a9ee --- /dev/null +++ b/python/remage/geombench/gdml_handling.py @@ -0,0 +1,260 @@ +# Copyright (C) 2025 Moritz Neuberger +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see . + + +from __future__ import annotations + +import tempfile +from pathlib import Path + +import pyg4ometry + + +def copy_referenced_resources( + lv: pyg4ometry.geant4.LogicalVolume, + src_reg: pyg4ometry.geant4.Registry, + dst_reg: pyg4ometry.geant4.Registry, +) -> None: + """Copy all resources referenced by a logical volume and its hierarchy. + + This function recursively copies materials, solids, positions, rotations, + and other defines needed by the logical volume and all its daughter volumes. + + Parameters + ---------- + lv : + The logical volume whose resources are to be copied. + src_reg : + The source registry from which to copy resources. + dst_reg : + The destination registry to which resources will be copied. + """ + + def _copy_named_object( + obj: object, + src_dict_name: str, + dst_dict_name: str, + src_reg: pyg4ometry.geant4.Registry, + dst_reg: pyg4ometry.geant4.Registry, + ) -> None: + """Copy a named object from source registry dict to destination if not present.""" + if not (obj and hasattr(obj, "name")): + return + + src_dict = getattr(src_reg, src_dict_name, {}) + if obj.name in src_dict and obj.name not in getattr(dst_reg, dst_dict_name, {}): + if not hasattr(dst_reg, dst_dict_name): + setattr(dst_reg, dst_dict_name, {}) + getattr(dst_reg, dst_dict_name)[obj.name] = src_dict[obj.name] + + def _copy_material( + material: pyg4ometry.geant4.Material, + src_reg: pyg4ometry.geant4.Registry, + dst_reg: pyg4ometry.geant4.Registry, + ) -> None: + """Copy a material and its property dependencies to the destination registry.""" + if not material or material.name in dst_reg.materialDict: + return + + # Copy material property references + if hasattr(material, "properties") and material.properties: + for prop_value in material.properties.values(): + if hasattr(prop_value, "name"): + for dict_name in ["matrixDict", "defineDict"]: + _copy_named_object( + prop_value, dict_name, dict_name, src_reg, dst_reg + ) + + material.registry = dst_reg + dst_reg.addMaterial(material) + + def _copy_solid_recursive( + solid: pyg4ometry.geant4.Solid, + src_reg: pyg4ometry.geant4.Registry, + dst_reg: pyg4ometry.geant4.Registry, + ) -> None: + """Recursively copy a solid and all its nested components.""" + if solid.name in dst_reg.solidDict: + return + + solid.registry = dst_reg + dst_reg.addSolid(solid) + + # Copy position lists (GenericPolycone, Tessellated, etc.) + for attr in ["pPositions", "vertices"]: + if hasattr(solid, attr) and getattr(solid, attr): + for pos in getattr(solid, attr): + _copy_named_object( + pos, "positionDict", "positionDict", src_reg, dst_reg + ) + + # Copy transforms for boolean solids + for tra_attr in ["tra1", "tra2"]: + if hasattr(solid, tra_attr) and getattr(solid, tra_attr): + tra = getattr(solid, tra_attr) + _copy_named_object( + tra.position, "positionDict", "positionDict", src_reg, dst_reg + ) + _copy_named_object( + tra.rotation, "rotationDict", "rotationDict", src_reg, dst_reg + ) + + # Recursively copy constituent solids for boolean operations + if hasattr(solid, "obj1") and solid.obj1: + _copy_solid_recursive(solid.obj1, src_reg, dst_reg) + if hasattr(solid, "obj2") and solid.obj2: + _copy_solid_recursive(solid.obj2, src_reg, dst_reg) + + _copy_material(lv.material, src_reg, dst_reg) + + if lv.solid: + _copy_solid_recursive(lv.solid, src_reg, dst_reg) + + # Copy define dictionary entries related to this volume + if hasattr(src_reg, "defineDict"): + for key in src_reg.defineDict: + if lv.name in key and key not in dst_reg.defineDict: + dst_reg.defineDict[key] = src_reg.defineDict[key] + + # Recursively copy daughter volumes + if hasattr(lv, "daughterVolumes"): + for daughter_pv in lv.daughterVolumes: + if hasattr(daughter_pv, "logicalVolume"): + daughter_lv = daughter_pv.logicalVolume + if daughter_lv.name not in dst_reg.logicalVolumeDict: + dst_reg.addLogicalVolume(daughter_lv) + copy_referenced_resources(daughter_lv, src_reg, dst_reg) + + +def load_gdml_geometry(gdml_path: Path, object_name: str = "object_lv") -> dict: + """Load a GDML geometry file and rename the world volume. + + Parameters + ---------- + gdml_path : + Path to the GDML file. + object_name : optional + Name to assign to the loaded world volume. Default is "object_lv". + + Returns + ------- + dict + Dictionary with keys: + - "object_lv": The loaded and renamed logical volume + - "registry": The pyg4ometry registry containing the geometry + """ + reader = pyg4ometry.gdml.Reader(gdml_path) + registry = reader.getRegistry() + world_volume = registry.getWorldVolume() + + # Rename the world volume + original_name = world_volume.name + world_volume.name = object_name + + # Update registry references + registry.logicalVolumeDict.pop(original_name, None) + registry.logicalVolumeList.remove(original_name) + registry.logicalVolumeDict[object_name] = world_volume + registry.logicalVolumeList.append(object_name) + + # Update world reference in registry + registry.setWorld(object_name) + + return {"object_lv": world_volume, "registry": registry} + + +def change_extent_of_world_volume( + geometry: dict, buffer_fraction: float = 0, object_name: str = "object_lv" +) -> dict: + """Expand the world volume to include buffer space around the geometry.""" + + world_lv = geometry["object_lv"] + registry = geometry["registry"] + + # Rename world volume if needed + if world_lv.name != object_name: + old_name = world_lv.name + world_lv.name = object_name + + # Update registry references + if old_name in registry.logicalVolumeDict: + registry.logicalVolumeDict.pop(old_name) + if old_name in registry.logicalVolumeList: + registry.logicalVolumeList.remove(old_name) + + registry.logicalVolumeDict[object_name] = world_lv + if object_name not in registry.logicalVolumeList: + registry.logicalVolumeList.append(object_name) + + # Always ensure world is set correctly + registry.setWorld(world_lv.name) + + # Calculate extent and create larger world box + extent = world_lv.extent(True) + width = (extent[1][0] - extent[0][0]) * (1 + buffer_fraction) + height = (extent[1][1] - extent[0][1]) * (1 + buffer_fraction) + depth = (extent[1][2] - extent[0][2]) * (1 + buffer_fraction) + + # Remove old solid from registry and create new larger box + old_solid_name = world_lv.solid.name + if old_solid_name in registry.solidDict: + del registry.solidDict[old_solid_name] + + new_box = pyg4ometry.geant4.solid.Box( + old_solid_name, width, height, depth, registry, "mm" + ) + + # Update the logical volume to use the new solid + world_lv.solid = new_box + + return {"object_lv": world_lv, "registry": registry} + + +def generate_tmp_gdml_geometry( + geometry: dict, buffer_fraction: float = 0.25, object_name: str = "object_lv" +) -> Path: + """Prepare a GDML geometry by wrapping it in a buffered world volume. + + Creates a new GDML file with the geometry positioned in a world volume + that includes buffer space around it. + + Parameters + ---------- + geometry : + Dictionary containing the loaded geometry with keys "object_lv" and "registry". + buffer_fraction : optional + Fractional buffer to add around the geometry. For example, 0.25 adds + 12.5%% extra space on each side. Default is 0.25. + object_name : optional + Name to assign to the object logical volume when positioning. + Default is "object_lv". + + Returns + ------- + Path + Path to the temporary GDML file with the adjusted geometry. + """ + positioned_geometry = change_extent_of_world_volume( + geometry, buffer_fraction=buffer_fraction, object_name=object_name + ) + + writer = pyg4ometry.gdml.Writer() + writer.addDetector(positioned_geometry["registry"]) + + temp_file = tempfile.NamedTemporaryFile(mode="w", delete=False, suffix=".gdml") # noqa: SIM115 + tempfile_path = Path(temp_file.name) + writer.write(tempfile_path) + + return tempfile_path diff --git a/python/remage/geombench/summary_generator.py b/python/remage/geombench/summary_generator.py new file mode 100644 index 00000000..80d1e83e --- /dev/null +++ b/python/remage/geombench/summary_generator.py @@ -0,0 +1,253 @@ +# Copyright (C) 2025 Moritz Neuberger +# +# This program is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the Free +# Software Foundation, either version 3 of the License, or (at your option) any +# later version. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +# details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with this program. If not, see . + + +from __future__ import annotations + +import argparse +from pathlib import Path + +import awkward as ak +import hist +import matplotlib.pyplot as plt +import numpy as np +import yaml +from lgdo import lh5 +from matplotlib.colors import LogNorm, Normalize + + +class SummaryGenerator: + """ + Class to generate summary analysis for remage geometry benchmark outputs. + """ + + def __init__( + self, sim_output_file: Path, args: argparse.Namespace, output_file_stem: str + ) -> None: + self.data_xy = lh5.read("benchmark_xy", sim_output_file).view_as("ak") + self.data_xz = lh5.read("benchmark_xz", sim_output_file).view_as("ak") + self.data_yz = lh5.read("benchmark_yz", sim_output_file).view_as("ak") + + self.n_x_gridpoint = len(np.unique(self.data_xy["x"])) + self.n_y_gridpoint = len(np.unique(self.data_xy["y"])) + self.n_z_gridpoint = len(np.unique(self.data_xz["z"])) + + self.x = np.linspace( + np.min(self.data_xy["x"]), + np.max(self.data_xy["x"]) + (self.data_xy["x"][1] - self.data_xy["x"][0]), + self.n_x_gridpoint + 1, + ) + self.y = np.linspace( + np.min(self.data_xy["y"]), + np.max(self.data_xy["y"]) + (self.data_yz["y"][1] - self.data_yz["y"][0]), + self.n_y_gridpoint + 1, + ) + self.z = np.linspace( + np.min(self.data_xz["z"]), + np.max(self.data_xz["z"]) + (self.data_xz["z"][1] - self.data_xz["z"][0]), + self.n_z_gridpoint + 1, + ) + + self.gdml_file = Path(args.geometry) + + self.output_file_stem = output_file_stem + self.output_dir = Path(args.output_dir) + + self.mult_map_3d = None + + def _multiplicative_reconstruction(self) -> None: + """ + Reconstruct 3D using multiplication of normalized projections. + This assumes independence and that hotspots must appear in all views. + """ + map_3d = np.ones((self.n_x_gridpoint, self.n_y_gridpoint, self.n_z_gridpoint)) + + # Normalize each projection + xy_norm = np.array(self.data_xy["time"]).reshape( + self.n_y_gridpoint, self.n_x_gridpoint + ) + xy_norm = xy_norm / np.max(xy_norm) if np.max(xy_norm) > 0 else xy_norm + + xz_norm = np.array(self.data_xz["time"]).reshape( + self.n_z_gridpoint, self.n_x_gridpoint + ) + xz_norm = xz_norm / np.max(xz_norm) if np.max(xz_norm) > 0 else xz_norm + + yz_norm = np.array(self.data_yz["time"]).reshape( + self.n_z_gridpoint, self.n_y_gridpoint + ) + yz_norm = yz_norm / np.max(yz_norm) if np.max(yz_norm) > 0 else yz_norm + + # Multiply projections (back-project and multiply) + for ix in range(self.n_x_gridpoint): + for iy in range(self.n_y_gridpoint): + for iz in range(self.n_z_gridpoint): + map_3d[ix, iy, iz] = ( + xy_norm[iy, ix] * xz_norm[iz, ix] * yz_norm[iz, iy] + ) + + self.mult_map_3d = map_3d + + def draw_simulation_time_profiles( + self, suffix: str = "simulation_time_profiles.pdf" + ) -> None: + """ + Draw the simulation times per event for the three projections. + """ + + def draw(self, norm): + fig, ax = plt.subplots(1, 3, figsize=(15, 5)) + # No self.fig; pass fig explicitly to helpers + + # Create histogram for XY plane + h_xy = hist.Hist( + hist.axis.Regular( + self.n_x_gridpoint, self.x[0], self.x[-1], name="x", label="x" + ), + hist.axis.Regular( + self.n_y_gridpoint, self.y[0], self.y[-1], name="y", label="y" + ), + ) + h_xy[:, :] = ( + np.array(self.data_xy["time"]) + .reshape(self.n_x_gridpoint, self.n_y_gridpoint) + .T + * 1e6 + ) + artists_xy = h_xy.plot2d( + ax=ax[0], norm=norm, edgecolor="face", linewidth=0, cbar=False + ) + ax[0].set_title("XY Plane") + ax[0].set_box_aspect((self.y[-1] - self.y[0]) / (self.x[-1] - self.x[0])) + cbar_xy = fig.colorbar(artists_xy[0], ax=ax[0], orientation="horizontal") + cbar_xy.set_label(r"Median sim time per event [$\mu$s]") + + # Create histogram for XZ plane + h_xz = hist.Hist( + hist.axis.Regular( + self.n_x_gridpoint, self.x[0], self.x[-1], name="x", label="x" + ), + hist.axis.Regular( + self.n_z_gridpoint, self.z[0], self.z[-1], name="z", label="z" + ), + ) + h_xz[:, :] = ( + np.array(self.data_xz["time"]) + .reshape(self.n_x_gridpoint, self.n_z_gridpoint) + .T + * 1e6 + ) + artists_xz = h_xz.plot2d( + ax=ax[1], norm=norm, edgecolor="face", linewidth=0, cbar=False + ) + ax[1].set_title("XZ Plane") + ax[1].set_box_aspect((self.z[-1] - self.z[0]) / (self.x[-1] - self.x[0])) + cbar_xz = fig.colorbar(artists_xz[0], ax=ax[1], orientation="horizontal") + cbar_xz.set_label(r"Median sim time per event [$\mu$s]") + + # Create histogram for YZ plane + h_yz = hist.Hist( + hist.axis.Regular( + self.n_y_gridpoint, self.y[0], self.y[-1], name="y", label="y" + ), + hist.axis.Regular( + self.n_z_gridpoint, self.z[0], self.z[-1], name="z", label="z" + ), + ) + h_yz[:, :] = ( + np.array(self.data_yz["time"]) + .reshape(self.n_y_gridpoint, self.n_z_gridpoint) + .T + * 1e6 + ) + artists_yz = h_yz.plot2d( + ax=ax[2], + norm=norm, + edgecolor="face", + linewidth=0, + cbar=False, + ) + ax[2].set_title("YZ Plane") + ax[2].set_box_aspect((self.z[-1] - self.z[0]) / (self.y[-1] - self.y[0])) + cbar_yz = fig.colorbar(artists_yz[0], ax=ax[2], orientation="horizontal") + cbar_yz.set_label(r"Median sim time per event [$\mu$s]") + + norm_string = "lin" + if isinstance(norm, LogNorm): + norm_string = "log" + + output_path = ( + self.output_dir / f"{self.output_file_stem}_{norm_string}_{suffix}" + ) + fig.savefig(output_path, bbox_inches="tight", dpi=300) + + plt.close(fig) + + draw(self, norm=LogNorm()) + draw(self, norm=Normalize()) + + def _get_hotspot_locations( + self, threshold: float = 0.8 + ) -> list[tuple[float, float, float]]: + """ + Get hotspot locations from multiplicative reconstruction above a certain threshold. + """ + if self.mult_map_3d is None: + self._multiplicative_reconstruction() + + hotspots = np.argwhere(self.mult_map_3d >= threshold * np.max(self.mult_map_3d)) + hotspot_coords = [] + for hs in hotspots: + x_coord = self.x[hs[0]] + y_coord = self.y[hs[1]] + z_coord = self.z[hs[2]] + hotspot_coords.append((float(x_coord), float(y_coord), float(z_coord))) + + return hotspot_coords + + def calculate_simulation_statistics(self, suffix: str = "_stats.yaml") -> dict: + """Calculate simulation statistics. + + Args: + suffix (str): Suffix for the output statistics file. + Returns: + dict: Dictionary containing simulation statistics. + + Details: + - Mean, standard deviation, minimum, and maximum simulation time per event. + - Hotspot locations identified from the multiplicative reconstruction. + + """ + + stats = { + "simulation_time_per_event": { + "mean": float(np.mean(ak.to_numpy(self.data_xy["time"]))), + "std": float(np.std(ak.to_numpy(self.data_xy["time"]))), + "min": float(np.min(ak.to_numpy(self.data_xy["time"]))), + "max": float(np.max(ak.to_numpy(self.data_xy["time"]))), + }, + "hotspots": self._get_hotspot_locations(threshold=0.8), + } + + output_path = self.output_dir / f"{self.output_file_stem}{suffix}" + + with Path(output_path).open("w") as f: + yaml.dump(stats, f) + + return stats + + def perform_analysis(self) -> dict: + self.draw_simulation_time_profiles() + return self.calculate_simulation_statistics() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 41ba6c56..77fcc249 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,10 +5,12 @@ set(_root ${PROJECT_SOURCE_DIR}) set(PROJECT_PUBLIC_HEADERS ${_root}/include/RMGAnalysisReader.hh + ${_root}/include/RMGGeomBenchOutputScheme.hh ${_root}/include/RMGDetectorHit.hh ${_root}/include/RMGDetectorMetadata.hh ${_root}/include/RMGExceptionHandler.hh ${_root}/include/RMGEventAction.hh + ${_root}/include/RMGGeomBench.hh ${_root}/include/RMGGeneratorCosmicMuons.hh ${_root}/include/RMGGeneratorFromFile.hh ${_root}/include/RMGGeneratorMUSUNCosmicMuons.hh @@ -55,11 +57,13 @@ set(PROJECT_PUBLIC_HEADERS set(PROJECT_SOURCES ${_root}/src/RMGAnalysisReader.cc + ${_root}/src/RMGGeomBenchOutputScheme.cc ${_root}/src/RMGDetectorHit.cc ${_root}/src/RMGExceptionHandler.cc ${_root}/src/RMGHardware.cc ${_root}/src/RMGHardwareMessenger.cc ${_root}/src/RMGEventAction.cc + ${_root}/src/RMGGeomBench.cc ${_root}/src/RMGGeneratorCosmicMuons.cc ${_root}/src/RMGGeneratorFromFile.cc ${_root}/src/RMGGeneratorMUSUNCosmicMuons.cc diff --git a/src/RMGGeomBench.cc b/src/RMGGeomBench.cc new file mode 100644 index 00000000..6006aa5c --- /dev/null +++ b/src/RMGGeomBench.cc @@ -0,0 +1,552 @@ +// Copyright (C) 2025 Moritz Neuberger +#include "RMGGeomBench.hh" + +#include +#include +#include + +#include "G4AnalysisManager.hh" +#include "G4GenericMessenger.hh" +#include "G4ParticleGun.hh" +#include "G4ParticleMomentum.hh" +#include "G4ParticleTypes.hh" +#include "G4ThreeVector.hh" +#include "G4TransportationManager.hh" +#include "G4VisExtent.hh" +#include "Randomize.hh" + +#include "RMGGeomBenchOutputScheme.hh" +#include "RMGHardware.hh" +#include "RMGLog.hh" +#include "RMGManager.hh" +#include "RMGNavigationTools.hh" + +namespace u = CLHEP; + +RMGGeomBench::RMGGeomBench() : RMGVGenerator("Benchmark") { + this->DefineCommands(); + fGun = std::make_unique(); + + // Set default values for configurable parameters + // negative means auto-calculate (30 pixels default for increment, world bounds for width) + user_increment = G4ThreeVector(-1.0, -1.0, -1.0); + sampling_width = G4ThreeVector(-1.0, -1.0, -1.0); +} + + +RMGGeomBench::~RMGGeomBench() = default; + +// Helper to find the benchmark output scheme if it's active +RMGGeomBenchOutputScheme* RMGGeomBench::GetBenchmarkOutputScheme() { + auto det_cons = RMGManager::Instance()->GetDetectorConstruction(); + for (auto& scheme : det_cons->GetAllActiveOutputSchemes()) { + auto benchmark_scheme = dynamic_cast(scheme.get()); + if (benchmark_scheme) return benchmark_scheme; + } + return nullptr; +} + +void RMGGeomBench::DefineCommands() { + + fMessenger = std::make_unique( + this, + "/RMG/Generator/Benchmark/", + "Commands for controlling the benchmarking simulation" + ); + + fMessenger->DeclarePropertyWithUnit("IncrementX", "mm", user_increment[0]) + .SetGuidance("Step size (increment) in X direction (negative = auto, default 30 pixels)") + .SetParameterName("dx", false) + .SetDefaultValue("-1.0"); + + fMessenger->DeclarePropertyWithUnit("IncrementY", "mm", user_increment[1]) + .SetGuidance("Step size (increment) in Y direction (negative = auto, default 30 pixels)") + .SetParameterName("dy", false) + .SetDefaultValue("-1.0"); + + fMessenger->DeclarePropertyWithUnit("IncrementZ", "mm", user_increment[2]) + .SetGuidance("Step size (increment) in Z direction (negative = auto, default 30 pixels)") + .SetParameterName("dz", false) + .SetDefaultValue("-1.0"); + + fMessenger->DeclarePropertyWithUnit("SamplingWidthX", "mm", sampling_width[0]) + .SetGuidance("Sampling width in X direction (negative = auto from world)") + .SetParameterName("wx", false) + .SetDefaultValue("-1.0"); + + fMessenger->DeclarePropertyWithUnit("SamplingWidthY", "mm", sampling_width[1]) + .SetGuidance("Sampling width in Y direction (negative = auto from world)") + .SetParameterName("wy", false) + .SetDefaultValue("-1.0"); + + fMessenger->DeclarePropertyWithUnit("SamplingWidthZ", "mm", sampling_width[2]) + .SetGuidance("Sampling width in Z direction (negative = auto from world)") + .SetParameterName("wz", false) + .SetDefaultValue("-1.0"); +} + + +void RMGGeomBench::BeginOfRunAction(const G4Run* r) { + + // Things set up at the begin of run action: + // - Start time of entire run (as usual) + // - Number of events to be processed + // - The NTuples for storing positional and timing information + // - The positional limits for the generator + // - Making sure the user doesn't break the program + + // In this current implementation, these parameters are hardcoded. + // Need to have a meeting with someone who understands the Python + // wrapper in order to get these into the CLI. + + // All lengths in mm + + // Initialize default geometry sampling parameters early to avoid + // using uninitialized member variables (was causing undefined behaviour). + cubesize = 10000.; + npixelsperrow = 30; // Legacy variable, kept for backward compatibility + + origin = G4ThreeVector(0., 0., 0.); + limit = origin - G4ThreeVector(cubesize / 2., cubesize / 2., cubesize / 2.); + + G4ThreeVector minn, maxx; // Only used for bounding conditions + + // Specify volume included in the benchmarking simulation (optional) + std::string targetvolumename = ""; + + // NOTE: The benchmark output scheme is registered as an optional output scheme + // in RMGUserInit::RegisterDefaultOptionalOutputSchemes() and must be activated + // by the user with: /RMG/Output/ActivateOutputScheme Benchmark + // It will be initialized by RMGRunAction::SetupAnalysisManager() if active. + + // whichntuple will be initialized to 0 on event 0 + // Start the timer + starttime = std::clock(); + currenttime = std::clock(); + + // Store total number of events to be processed + totalnevents = r->GetNumberOfEventToBeProcessed(); + + // The benchmarker assumes the world volume is a cube when calculating bounds + auto worldsolid = G4TransportationManager::GetTransportationManager() + ->GetNavigatorForTracking() + ->GetWorldVolume() + ->GetLogicalVolume() + ->GetSolid(); + auto typeofworld = worldsolid->GetEntityType(); + if (typeofworld != "G4Box") + RMGLog::OutFormat( + RMGLog::warning, + "World entity shape ({}) is not a G4Box (why not? :| ). A very large benchmarking cube may " + "crash the simulation.", + typeofworld + ); + + // set the origin and limit according to the world volume + G4ThreeVector min, max; + worldsolid->BoundingLimits(min, max); + + // Set up sampling dimensions based on user configuration or world bounds + if (sampling_width.x() < 0) { sampling_width.setX(max.x() - min.x()); } + if (sampling_width.y() < 0) { sampling_width.setY(max.y() - min.y()); } + if (sampling_width.z() < 0) { sampling_width.setZ(max.z() - min.z()); } + + // Calculate increments and number of pixels + // If user specified increment, use it; otherwise default to 30 pixels + if (user_increment.x() > 0) { + increment.setX(user_increment.x()); + npixels_x = static_cast(std::ceil(sampling_width.x() / increment.x())); + } else { + npixels_x = 30; // default + increment.setX(sampling_width.x() / npixels_x); + } + + if (user_increment.y() > 0) { + increment.setY(user_increment.y()); + npixels_y = static_cast(std::ceil(sampling_width.y() / increment.y())); + } else { + npixels_y = 30; // default + increment.setY(sampling_width.y() / npixels_y); + } + + if (user_increment.z() > 0) { + increment.setZ(user_increment.z()); + npixels_z = static_cast(std::ceil(sampling_width.z() / increment.z())); + } else { + npixels_z = 30; // default + increment.setZ(sampling_width.z() / npixels_z); + } + + // Center the sampling region on the world volume center + origin = (max + min) / 2.; + limit = origin - sampling_width / 2.; + + + RMGLog::Out(RMGLog::summary, "Benchmark sampling configuration:"); + RMGLog::OutFormat( + RMGLog::summary, + " X: {} pixels, width = {} mm, increment = {} mm", + npixels_x, + sampling_width.x(), + increment.x() + ); + RMGLog::OutFormat( + RMGLog::summary, + " Y: {} pixels, width = {} mm, increment = {} mm", + npixels_y, + sampling_width.y(), + increment.y() + ); + RMGLog::OutFormat( + RMGLog::summary, + " Z: {} pixels, width = {} mm, increment = {} mm", + npixels_z, + sampling_width.z(), + increment.z() + ); + RMGLog::OutFormat( + RMGLog::summary, + "Sampling region origin at (X,Y,Z): ({}, {}, {}) mm", + origin.x(), + origin.y(), + origin.z() + ); + RMGLog::OutFormat( + RMGLog::summary, + "Sampling region lower limit at (X,Y,Z): ({}, {}, {}) mm", + limit.x(), + limit.y(), + limit.z() + ); + + // Calculate the number of pixels to simulate, and number of primaries per pixel + // NOW that we know npixels_x, npixels_y, npixels_z + + // Total pixels across all three planes (XZ, YZ, XY) + int pixels_xz = npixels_x * npixels_z; + int pixels_yz = npixels_y * npixels_z; + int pixels_xy = npixels_x * npixels_y; + totalnpixels = pixels_xz + pixels_yz + pixels_xy; + + if (totalnpixels <= 0) { + RMGLog::Out(RMGLog::fatal, "Total number of pixels is zero or negative! Exiting..."); + } + + neventsperpixel = totalnevents / totalnpixels; + + // Calculate events per bunch (1% of events per pixel) + events_per_bunch = std::max(1, static_cast(std::ceil(neventsperpixel * 0.01))); + + // Calculate total number of batch rounds + total_batch_rounds = static_cast( + std::ceil(static_cast(neventsperpixel) / events_per_bunch) + ); + + RMGLog::OutFormat( + RMGLog::summary, + "Total pixels to sample: {} (XZ: {}, YZ: {}, XY: {})", + totalnpixels, + pixels_xz, + pixels_yz, + pixels_xy + ); + RMGLog::OutFormat(RMGLog::debug, "Events per bunch: {}", events_per_bunch); + RMGLog::OutFormat(RMGLog::debug, "Total batch rounds: {}", total_batch_rounds); + + // Safeties + + // Could maybe downgrade this from a warning to a summary... + if (neventsperpixel != floor(neventsperpixel)) { + RMGLog::OutFormat( + RMGLog::warning, + "Specified number of primaries({}) doesn't divide evenly into the specified number of " + "pixels ({})", + totalnevents, + totalnpixels + ); + RMGLog::Out(RMGLog::warning, "Rounding down to nearest integer number of primaries per pixel."); + } + neventsperpixel = floor(neventsperpixel); + + // Having a single event sampled per pixel breaks the generator algorithm + // Doing this would also be completely useless from a benchmarking perspective + if (neventsperpixel < 2) + RMGLog::Out(RMGLog::fatal, "Not enough primaries to sample each pixel at least twice! Exiting..."); + else RMGLog::OutFormat(RMGLog::summary, "Number of primaries per pixel: {}", neventsperpixel); + +} // BeginOfRunAction + + +void RMGGeomBench::EndOfRunAction(const G4Run* /*r*/) { + // Save the last batch if there's any data + if (current_batch_event > 0) { + double bunch_time = static_cast(std::clock() - bunchstarttime) / CLOCKS_PER_SEC; + RecordBatchTime(current_pixel_index, bunch_time); + RMGLog::Out(RMGLog::debug, "Recorded final batch time in EndOfRunAction"); + } + + // Save all pixel data + RMGLog::Out(RMGLog::debug, "Saving all pixel data in EndOfRunAction"); + SaveAllPixels(); +} // EndOfRunAction + + +void RMGGeomBench::RecordBatchTime(size_t pixel_idx, double batch_time) { + if (pixel_idx >= totalnpixels) { + RMGLog::OutFormat(RMGLog::warning, "Invalid pixel index ({}) in RecordBatchTime", pixel_idx); + return; + } + + // Ensure the vector is large enough + if (static_cast(pixel_batch_times.size()) <= pixel_idx) { + pixel_batch_times.resize(pixel_idx + 1); + } + + pixel_batch_times[pixel_idx].push_back(batch_time); + + RMGLog::OutFormat(RMGLog::debug, "Recorded batch time {} s for pixel {}", batch_time, pixel_idx); +} // RecordBatchTime + +void RMGGeomBench::SaveAllPixels() { + // Find the benchmark output scheme if it's active + auto benchmark_scheme = GetBenchmarkOutputScheme(); + if (!benchmark_scheme) { + RMGLog::OutDev( + RMGLog::warning, + "Benchmark output scheme not active. Use /RMG/Output/ActivateOutputScheme Benchmark" + ); + return; + } + + // Iterate through all pixels and save their data + int pixel_idx = 0; + + for (int plane = 0; plane < 3; plane++) { + int max_i = 0, max_j = 0; + + switch (plane) { + case 0: // XZ plane + max_i = npixels_x; + max_j = npixels_z; + break; + case 1: // YZ plane + max_i = npixels_y; + max_j = npixels_z; + break; + case 2: // XY plane + max_i = npixels_x; + max_j = npixels_y; + break; + } + + for (int j = 0; j < max_j; j++) { + for (int i = 0; i < max_i; i++) { + // Calculate position for this pixel + double x_pos, y_pos, z_pos; + + switch (plane) { + case 0: // XZ plane + x_pos = limit.x() + i * increment.x(); + y_pos = limit.y(); + z_pos = limit.z() + j * increment.z(); + break; + case 1: // YZ plane + x_pos = limit.x(); + y_pos = limit.y() + i * increment.y(); + z_pos = limit.z() + j * increment.z(); + break; + case 2: // XY plane + x_pos = limit.x() + i * increment.x(); + y_pos = limit.y() + j * increment.y(); + z_pos = limit.z(); + break; + } + + // Calculate median time from batch times + double median_time_per_event = 0.0; + + if (pixel_idx < static_cast(pixel_batch_times.size()) && + !pixel_batch_times[pixel_idx].empty()) { + std::vector sorted_times = pixel_batch_times[pixel_idx]; + std::sort(sorted_times.begin(), sorted_times.end()); + + size_t n = sorted_times.size(); + double median_time; + if (n % 2 == 0) { + median_time = (sorted_times[n / 2 - 1] + sorted_times[n / 2]) / 2.0; + } else { + median_time = sorted_times[n / 2]; + } + + if (events_per_bunch > 0) { median_time_per_event = median_time / events_per_bunch; } + + RMGLog::OutFormat( + RMGLog::debug, + "Pixel {} (plane {}): {} batches, median time per event: {} s", + pixel_idx, + plane, + n, + median_time_per_event + ); + } + + benchmark_scheme->SavePixel(plane, x_pos, y_pos, z_pos, median_time_per_event); + pixel_idx++; + } + } + } + + RMGLog::OutFormat(RMGLog::summary, "Saved data for all {} pixels", pixel_idx); +} // SaveAllPixels + + +void RMGGeomBench::GeneratePrimaries(G4Event* event) { + + ID = event->GetEventID(); + + if (ID >= neventsperpixel * totalnpixels) + return; // Should only apply in situations where the nevents doesn't divide evenly into the npixels + + if (ID == 0) { + // Originally initialized in BeginOfRunAction, but delay between that and evt 0 is non-trivial + currenttime = std::clock(); + bunchstarttime = std::clock(); + + current_batch_event = 0; + current_pixel_index = 0; + current_batch_round = 0; + + // Initialize storage for batch times + pixel_batch_times.clear(); + pixel_batch_times.resize(totalnpixels); + } + + // Calculate which batch round and which pixel we're in + // Structure: batch_round -> pixel -> events within batch + // Total events = total_batch_rounds * totalnpixels * events_per_bunch (approximately) + int events_so_far = ID; + current_batch_round = events_so_far / (totalnpixels * events_per_bunch); + int remainder = events_so_far % (totalnpixels * events_per_bunch); + current_pixel_index = remainder / events_per_bunch; + current_batch_event = remainder % events_per_bunch; + + // Check if we're starting a new batch + if (ID == 0 || current_batch_event == 0) { + // Save the previous batch time if this isn't the first event + if (ID > 0) { + double bunch_time = static_cast(std::clock() - bunchstarttime) / CLOCKS_PER_SEC; + int prev_pixel_index = (ID - 1) % (totalnpixels * events_per_bunch) / events_per_bunch; + RecordBatchTime(prev_pixel_index, bunch_time); + RMGLog::OutFormat( + RMGLog::debug, + "Batch complete for pixel {} at event # {}, CPU time: {} s", + prev_pixel_index, + ID - 1, + bunch_time + ); + } + + // Start timing for new batch + bunchstarttime = std::clock(); + + if (current_batch_event == 0 && current_pixel_index == 0) { + RMGLog::OutFormat(RMGLog::debug, "Starting batch round {} at event # {}", current_batch_round, ID); + } + } + + // Convert linear pixel index to plane and coordinates + int pixels_xz = npixels_x * npixels_z; + int pixels_yz = npixels_y * npixels_z; + + int plane; + int i_index, j_index; + + if (current_pixel_index < pixels_xz) { + // XZ plane + plane = 0; + i_index = current_pixel_index % npixels_x; + j_index = current_pixel_index / npixels_x; + } else if (current_pixel_index < pixels_xz + pixels_yz) { + // YZ plane + plane = 1; + int local_idx = current_pixel_index - pixels_xz; + i_index = local_idx % npixels_y; + j_index = local_idx / npixels_y; + } else { + // XY plane + plane = 2; + int local_idx = current_pixel_index - pixels_xz - pixels_yz; + i_index = local_idx % npixels_x; + j_index = local_idx / npixels_x; + } + + // Calculate position for this pixel + G4ThreeVector current_position; + switch (plane) { + case 0: // XZ plane + current_position.setX(limit.x() + i_index * increment.x()); + current_position.setY(limit.y()); + current_position.setZ(limit.z() + j_index * increment.z()); + break; + case 1: // YZ plane + current_position.setX(limit.x()); + current_position.setY(limit.y() + i_index * increment.y()); + current_position.setZ(limit.z() + j_index * increment.z()); + break; + case 2: // XY plane + current_position.setX(limit.x() + i_index * increment.x()); + current_position.setY(limit.y() + j_index * increment.y()); + current_position.setZ(limit.z()); + break; + } + + fGun->SetNumberOfParticles(1); + fGun->SetParticleDefinition(G4Geantino::Definition()); + + // Randomize particle position within the current grid space + // For each face of the sampling cube, we will have two + // randomized positional values, with the current position + // as a lower bound, and one fixed positional value, which + // is the face we're currently sampling from + // To make the code prettier, I'll just sample each + // value uniformly, then squish whichever one is constant + + double xtemp = G4UniformRand() * increment.x() + current_position.x(); + double ytemp = G4UniformRand() * increment.y() + current_position.y(); + double ztemp = G4UniformRand() * increment.z() + current_position.z(); + + G4ThreeVector momentumdir(0., 0., 0.); + + switch (plane) { + case 0: { // XZ + ytemp = limit.y(); + momentumdir.setY(1.); + break; + } + case 1: { // YZ + xtemp = limit.x(); + momentumdir.setX(1.); + break; + } + case 2: { // XY + ztemp = limit.z(); + momentumdir.setZ(1.); + break; + } + default: { + RMGLog::OutFormat(RMGLog::fatal, "Invalid plane ({}) in GeneratePrimaries()", plane); + } + } // switch(plane) + + fGun->SetParticlePosition({xtemp, ytemp, ztemp}); + fGun->SetParticleMomentumDirection(momentumdir); + fGun->SetParticleEnergy(1 * u::GeV); + + RMGLog::OutFormat(RMGLog::debug, "Sampled X: {} Sampled Y: {} Sampled Z: {}", xtemp, ytemp, ztemp); + RMGLog::Out(RMGLog::debug, "Momentum direction: ", momentumdir); + + fGun->GeneratePrimaryVertex(event); + +} // GeneratePrimaries + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/src/RMGGeomBenchOutputScheme.cc b/src/RMGGeomBenchOutputScheme.cc new file mode 100644 index 00000000..db490566 --- /dev/null +++ b/src/RMGGeomBenchOutputScheme.cc @@ -0,0 +1,105 @@ +// Copyright (C) 2025 Moritz Neuberger +// +// This program is free software: you can redistribute it and/or modify it under +// the terms of the GNU Lesser General Public License as published by the Free +// Software Foundation, either version 3 of the License, or (at your option) any +// later version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more +// details. +// +// You should have received a copy of the GNU Lesser General Public License +// along with this program. If not, see . + +#include "RMGGeomBenchOutputScheme.hh" + +#include "G4AnalysisManager.hh" + +#include "RMGLog.hh" +#include "RMGOutputManager.hh" + +RMGGeomBenchOutputScheme::RMGGeomBenchOutputScheme() { this->DefineCommands(); } + +// invoked in RMGRunAction::SetupAnalysisManager() +void RMGGeomBenchOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) { + + auto rmg_man = RMGOutputManager::Instance(); + + // Create auxiliary ntuple for XZ plane benchmark + fNtupleIDs[0] = rmg_man->CreateAndRegisterAuxNtuple("benchmark_xz", "RMGGeomBenchOutputScheme", ana_man); + ana_man->CreateNtupleDColumn(fNtupleIDs[0], "x"); + ana_man->CreateNtupleDColumn(fNtupleIDs[0], "z"); + ana_man->CreateNtupleDColumn(fNtupleIDs[0], "time"); + ana_man->FinishNtuple(fNtupleIDs[0]); + + // Create auxiliary ntuple for YZ plane benchmark + fNtupleIDs[1] = rmg_man->CreateAndRegisterAuxNtuple("benchmark_yz", "RMGGeomBenchOutputScheme", ana_man); + ana_man->CreateNtupleDColumn(fNtupleIDs[1], "y"); + ana_man->CreateNtupleDColumn(fNtupleIDs[1], "z"); + ana_man->CreateNtupleDColumn(fNtupleIDs[1], "time"); + ana_man->FinishNtuple(fNtupleIDs[1]); + + // Create auxiliary ntuple for XY plane benchmark + fNtupleIDs[2] = rmg_man->CreateAndRegisterAuxNtuple("benchmark_xy", "RMGGeomBenchOutputScheme", ana_man); + ana_man->CreateNtupleDColumn(fNtupleIDs[2], "x"); + ana_man->CreateNtupleDColumn(fNtupleIDs[2], "y"); + ana_man->CreateNtupleDColumn(fNtupleIDs[2], "time"); + ana_man->FinishNtuple(fNtupleIDs[2]); +} + +void RMGGeomBenchOutputScheme::SavePixel(int plane_id, double x, double y, double z, double time) { + + if (plane_id < 0 || plane_id > 2) { + RMGLog::Out(RMGLog::warning, "Invalid plane_id (", plane_id, ") in SavePixel(); skipping"); + return; + } + + auto rmg_man = RMGOutputManager::Instance(); + if (!rmg_man->IsPersistencyEnabled()) return; + + const auto ana_man = G4AnalysisManager::Instance(); + if (!ana_man) { + RMGLog::Out(RMGLog::warning, "Analysis manager not available in SavePixel(); skipping"); + return; + } + + int ntuple_id = fNtupleIDs[plane_id]; + if (ntuple_id < 0) { + RMGLog::Out(RMGLog::warning, "Ntuple ID not initialized for plane_id=", plane_id); + return; + } + + int col_id = 0; + switch (plane_id) { + case 0: // XZ plane + ana_man->FillNtupleDColumn(ntuple_id, col_id++, x); + ana_man->FillNtupleDColumn(ntuple_id, col_id++, z); + break; + case 1: // YZ plane + ana_man->FillNtupleDColumn(ntuple_id, col_id++, y); + ana_man->FillNtupleDColumn(ntuple_id, col_id++, z); + break; + case 2: // XY plane + ana_man->FillNtupleDColumn(ntuple_id, col_id++, x); + ana_man->FillNtupleDColumn(ntuple_id, col_id++, y); + break; + } + + ana_man->FillNtupleDColumn(ntuple_id, col_id++, time); + ana_man->AddNtupleRow(ntuple_id); +} + +void RMGGeomBenchOutputScheme::DefineCommands() { + + fMessenger = std::make_unique( + this, + "/RMG/Output/Benchmark/", + "Commands for controlling geometry navigation benchmark output." + ); + + // Future commands for configuring benchmark output could go here +} + +// vim: tabstop=2 shiftwidth=2 expandtab diff --git a/src/RMGMasterGenerator.cc b/src/RMGMasterGenerator.cc index e042b130..08a7c7e5 100644 --- a/src/RMGMasterGenerator.cc +++ b/src/RMGMasterGenerator.cc @@ -27,6 +27,7 @@ #include "RMGGeneratorG4Gun.hh" #include "RMGGeneratorGPS.hh" #include "RMGGeneratorMUSUNCosmicMuons.hh" +#include "RMGGeomBench.hh" #include "RMGLog.hh" #include "RMGManager.hh" #include "RMGTools.hh" @@ -124,6 +125,7 @@ void RMGMasterGenerator::SetGenerator(RMGMasterGenerator::Generator gen) { fGeneratorObj = std::make_unique(); break; case Generator::kFromFile: fGeneratorObj = std::make_unique(); break; + case Generator::kGeomBench: fGeneratorObj = std::make_unique(); break; case Generator::kUndefined: case Generator::kUserDefined: break; default: diff --git a/src/remage-doc-dump.cc b/src/remage-doc-dump.cc index 515e69cb..34f7eba4 100644 --- a/src/remage-doc-dump.cc +++ b/src/remage-doc-dump.cc @@ -30,6 +30,7 @@ #endif #include "RMGGeneratorFromFile.hh" #include "RMGGeneratorMUSUNCosmicMuons.hh" +#include "RMGGeomBenchOutputScheme.hh" #include "RMGGermaniumOutputScheme.hh" #include "RMGInnerBremsstrahlungProcess.hh" #include "RMGIsotopeFilterScheme.hh" @@ -63,6 +64,7 @@ void init_extra() { new RMGIsotopeFilterScheme(); new RMGTrackOutputScheme(); new RMGParticleFilterScheme(); + new RMGGeomBenchOutputScheme(); // confinments new RMGVertexConfinement(); diff --git a/tests/python/test_geombench_cli.py b/tests/python/test_geombench_cli.py new file mode 100644 index 00000000..f21e8516 --- /dev/null +++ b/tests/python/test_geombench_cli.py @@ -0,0 +1,56 @@ +from __future__ import annotations + +import re +from pathlib import Path + +from lgdo import lh5 +from remage.geombench.cli import generate_macro, remage_geombench_cli + + +def test_remage_geombench_cli(tmp_path: Path) -> None: + geometry_path = "gdml/geometry.gdml" + + output_dir = tmp_path / "output" + output_dir.mkdir() + + args = [ + geometry_path, + "--output-dir", + str(output_dir), + "--num-events", + "100000", + "--grid-increment", + "1", + ] + + remage_geombench_cli(external_args=args) + + assert (output_dir / "geometry.lh5").exists() + assert (output_dir / "geometry_lin_simulation_time_profiles.pdf").exists() + assert (output_dir / "geometry_stats.yaml").exists() + + # check if data was written correctly + assert lh5.read("benchmark_xy", output_dir / "geometry.lh5") is not None + + +def test_generate_macro(tmp_path: Path) -> None: + class Args: + geometry = Path("gdml/geometry.gdml") + output_dir = tmp_path / "output" + num_events = 1000000 + grid_increment = 2 + grid_increments = "" + dry_run = False + + args = Args() + + output_file_stem = args.geometry.stem + + macro_content = generate_macro(args, output_file_stem=output_file_stem) + + assert re.search( + r"/RMG/Geometry/IncludeGDMLFile .*/gdml/geometry.gdml", macro_content + ) + assert re.search(r"/RMG/Output/FileName .*/output/geometry.lh5", macro_content) + assert re.search(r"/RMG/Generator/Benchmark/IncrementX 2", macro_content) + assert re.search(r"/run/beamOn 1000000", macro_content) diff --git a/tests/python/test_geombench_gdml.py b/tests/python/test_geombench_gdml.py new file mode 100644 index 00000000..88296139 --- /dev/null +++ b/tests/python/test_geombench_gdml.py @@ -0,0 +1,22 @@ +from __future__ import annotations + +from pathlib import Path + +from remage.geombench.gdml_handling import ( + generate_tmp_gdml_geometry, + load_gdml_geometry, +) + + +def test_generate_tmp_gdml_geometry() -> None: + geometry_path = "gdml/geometry.gdml" + component_name = "germanium" + + loaded_gdml = load_gdml_geometry(geometry_path) + + generated_file_path = generate_tmp_gdml_geometry(loaded_gdml) + + with Path(generated_file_path).open("r") as f: + gdml_content = f.read() + + assert f'name="{component_name}"' in gdml_content