Initial commit

CMakeLists.txt

+cmake_minimum_required(VERSION 3.12 FATAL_ERROR)
+cmake_policy(VERSION 3.12...3.18)
+
+project(g4run)
+
+find_package(Geant4 REQUIRED gdml)
+
+add_executable(g4run
+  actions.cc
+  field.cc
+  geometry.cc
+  init.cc
+  main.cc
+  primary.cc
+)
+
+separate_arguments(Geant4_CXX_FLAGS)
+
+target_compile_options(g4run PRIVATE ${Geant4_CXX_FLAGS})
+target_compile_definitions(g4run PRIVATE ${Geant4_DEFINITIONS})
+target_include_directories(g4run PRIVATE ${Geant4_INCLUDE_DIRS})
+target_link_libraries(g4run PRIVATE ${Geant4_LIBRARIES})
+
+include(GNUInstallDirs)
+
+install(TARGETS g4run RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})

LICENSE

+BSD 2-Clause License
+
+Copyright (c) 2020, Guilherme Amadio
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

README.md

+# G4 - Simple Generic Program to Run Geant4 Simulations
+
+This program takes as input a GDML file, a reference physics list,
+a primary particle type and its initial energy, and runs an event
+loop in batch mode or in an interactive session. Its intented use
+is in profiling and benchmarking Geant4 performance.

actions.cc

+#include "actions.h"
+
+static bool stats = false;
+
+void enable_statistics() { stats = true; }
+
+G4Run* RunAction::GenerateRun() { return nullptr; }
+void RunAction::BeginOfRunAction(const G4Run*) { }
+void RunAction::EndOfRunAction(const G4Run*) { }
+
+void EventAction::BeginOfEventAction(const G4Event*) { }
+void EventAction::EndOfEventAction(const G4Event*) { }
+
+void TrackingAction::PreUserTrackingAction(const G4Track*) { }
+void TrackingAction::PostUserTrackingAction(const G4Track*) { }
+
+void SteppingAction::UserSteppingAction(const G4Step*) { }

actions.h

+#ifndef ACTIONS_H
+#define ACTIONS_H
+
+#include <G4UserRunAction.hh>
+#include <G4UserEventAction.hh>
+#include <G4UserTrackingAction.hh>
+#include <G4UserSteppingAction.hh>
+
+class G4Run;
+class G4Event;
+class G4Track;
+class G4Step;
+
+void enable_statistics();
+
+class RunAction final : public G4UserRunAction {
+public:
+	G4Run* GenerateRun() override;
+	void BeginOfRunAction(const G4Run*) override;
+	void EndOfRunAction(const G4Run*) override;
+};
+
+class EventAction final : public G4UserEventAction {
+public:
+	void BeginOfEventAction(const G4Event*) override;
+	void EndOfEventAction(const G4Event*) override;
+};
+
+class SteppingAction final : public G4UserSteppingAction {
+public:
+	void UserSteppingAction(const G4Step*) override;
+};
+
+class TrackingAction final : public G4UserTrackingAction {
+public:
+	void PreUserTrackingAction(const G4Track*) override;
+	void PostUserTrackingAction(const G4Track*) override;
+};
+
+#endif

field.cc

+#include <G4FieldManager.hh>
+#include <G4SystemOfUnits.hh>
+#include <G4TransportationManager.hh>
+#include <G4UniformMagField.hh>
+
+static inline G4FieldManager *field_manager()
+{
+	return G4TransportationManager::GetTransportationManager()->GetFieldManager();
+}
+
+void set_uniform_magnetic_field(double Bx, double By, double Bz)
+{
+	auto magnetic_field = G4ThreeVector(Bx, By, Bz) * tesla;
+
+	if (magnetic_field.mag2() > 0.0) {
+		auto field = new G4UniformMagField(magnetic_field);
+		field_manager()->SetDetectorField(field);
+		field_manager()->CreateChordFinder(field);
+	}
+}

field.h

+#ifndef _FIELD_H
+#define _FIELD_H
+
+void set_uniform_magnetic_field(double Bx, double By, double Bz);
+
+#endif

geometry.cc

+#include "geometry.h"
+
+#include <G4GDMLParser.hh>
+
+#include <err.h>
+#include <unistd.h>
+
+G4VPhysicalVolume *DetectorConstruction::Construct()
+{
+	if (access(gdml, R_OK))
+		err(errno, "%s", gdml);
+
+	G4GDMLParser p;
+
+	p.Read(gdml, /* validate */ true);
+
+	auto world = p.GetWorldVolume();
+
+    if (!world)
+        errx(1, "error reading geometry file");
+
+    return world;
+}

geometry.h

+#ifndef _GEOMETRY_H
+#define _GEOMETRY_H
+
+#include <G4VUserDetectorConstruction.hh>
+
+class DetectorConstruction final : public G4VUserDetectorConstruction {
+public:
+    DetectorConstruction() = delete;
+    DetectorConstruction(const char* path) : gdml(path) {}
+	G4VPhysicalVolume *Construct() override;
+private:
+    const char *gdml;
+};
+
+#endif

init.cc

+#include "init.h"
+
+#include "actions.h"
+#include "primary.h"
+
+void InitializationAction::Build() const
+{
+  SetUserAction(new RunAction());
+  SetUserAction(new PrimaryGeneratorAction());
+}

init.h

+#ifndef _INIT_H
+#define _INIT_H
+
+#include <G4VUserActionInitialization.hh>
+
+class InitializationAction final : public G4VUserActionInitialization {
+public:
+	void Build() const override;
+};
+
+#endif

main.cc

+#include <G4Version.hh>
+
+#if !defined(G4MULTITHREADED)
+#include <G4RunManager.hh>
+using RunManager = G4RunManager;
+#else
+#include <G4MTRunManager.hh>
+using RunManager = G4MTRunManager;
+#endif
+
+#include <G4UImanager.hh>
+#include <G4PhysListFactory.hh>
+
+#include "init.h"
+#include "field.h"
+#include "actions.h"
+#include "primary.h"
+#include "geometry.h"
+
+#include <cstdio>
+
+#include <err.h>
+#include <libgen.h>
+#include <getopt.h>
+#include <unistd.h>
+
+/* command line options */
+
+struct option options[] = {
+	{ "help",         no_argument,       NULL, 'h'},
+	{ "interactive",  no_argument,       NULL, 'i'},
+	{ "quiet",        no_argument,       NULL, 'q'},
+	{ "stats",        no_argument,       NULL, 'S'},
+	{ "verbose",      no_argument,       NULL, 'v'},
+	{ "version",      no_argument,       NULL, 'V'},
+	{ "events",       required_argument, NULL, 'e'},
+	{ "energy",       required_argument, NULL, 'E'},
+	{ "field",        required_argument, NULL, 'f'},
+	{ "geometry",     required_argument, NULL, 'g'},
+	{ "macro",        required_argument, NULL, 'm'},
+	{ "particle",     required_argument, NULL, 'p'},
+	{ "physics-list", required_argument, NULL, 'l'},
+	{ "seed",         required_argument, NULL, 's'},
+	{ "threads",      required_argument, NULL, 't'},
+};
+
+/* globals */
+
+int nthreads = 1;
+int verbose = false;
+bool interactive = false;
+long long seed = 0;
+long long nevents = 0;
+const char *physics_list = "FTFP_BERT";
+char *gdml_filename = nullptr;
+char *macro = nullptr;
+double Bz = 0.0; /* in Tesla */
+
+void help(const char *name)
+{
+	printf("Usage: %s [OPTIONS]\n\n", name);
+
+	for (auto& opt : options)
+		printf("\t-%c  --%-20s%-20s\n\n",
+			opt.val, opt.name, opt.has_arg ? opt.name : "");
+}
+
+void parse_options(int argc, char **argv)
+{
+	if (argc == 1) {
+		help(basename(argv[0]));
+		exit(0);
+	}
+
+	int optidx = 0;
+
+	for(;;) {
+		switch(getopt_long(argc, argv, "hiqSvVe:E:f:g:m:p:l:s:t:", options, &optidx)) {
+		case -1:
+			return;
+
+		case 'h':
+			help(basename(argv[0]));
+			exit(0);
+
+		case 'i':
+			interactive = true;
+			break;
+
+		case 'q':
+			verbose = 0;
+			break;
+
+		case 'v':
+			++verbose;
+			break;
+
+		case 'V': {
+			const char *start = strchr(G4Version, '-');
+			const char *stop  = strrchr(G4Version, '$');
+			char version[64] = { 0, };
+			strncpy(version, start+1, (size_t)(stop - start - 1));
+			printf("%s for Geant4 %s\n", basename(argv[0]), version);
+			exit(0);
+		}
+
+		case 'e':
+			nevents = (long long) strtol(optarg, nullptr, 10);
+
+			if (nevents < 0)
+				errx(EINVAL, "invalid number of events: %lld", nevents);
+			break;
+
+		case 'E': {
+			double E = (double) strtod(optarg, nullptr);
+			if (E >= 0.0)
+				set_primary_energy(E);
+			else
+				errx(EINVAL, "invalid value for the energy: %s", optarg);
+			break;
+		}
+
+		case 'f':
+			Bz = (double) strtod(optarg, nullptr);
+			break;
+
+		case 'g':
+			gdml_filename = optarg;
+			break;
+
+		case 'm':
+			macro = optarg;
+			break;
+
+		case 'p':
+			set_primary_name(optarg);
+			break;
+
+		case 'l':
+			physics_list = optarg;
+			break;
+
+		case 's':
+			seed = strtoll(optarg, nullptr, 10);
+			break;
+
+		case 'S':
+			enable_statistics();
+			break;
+
+		case 't':
+#ifdef G4MULTITHREADED
+			nthreads = (int) strtol(optarg, nullptr, 10);
+
+			if (nthreads <= 0)
+				errx(EINVAL, "invalid number of threads");
+#else
+			warnx("Geant4 has multithreading disabled");
+#endif
+			break;
+
+		default:
+			help(basename(argv[0]));
+			errx(1, "unknown option");
+		}
+	}
+}
+
+static const char *tty;
+
+void disable_stdout()
+{
+	tty = ttyname(fileno(stdout));
+	if(!freopen("/dev/null", "a+", stdout))
+		errx(1, "failed to disable stdout");
+}
+
+void enable_stdout()
+{
+	if(!freopen(tty, "w", stdout))
+		errx(1, "failed to enable stdout");
+}
+
+int main(int argc, char **argv)
+{
+	parse_options(argc, argv);
+
+	if (!verbose)
+		disable_stdout();
+
+	if (seed)
+		G4Random::setTheSeed(seed);
+
+	if (!gdml_filename)
+		errx(1, "no geometry file specified");
+
+	G4PhysListFactory factory;
+
+	if (!factory.IsReferencePhysList(physics_list))
+		errx(1, "unknown physics list: %s", physics_list);
+
+	auto runManager = new RunManager();
+
+	G4UImanager* UI = G4UImanager::GetUIpointer();
+
+	UI->ApplyCommand(G4String("/process/verbose       0"));
+	UI->ApplyCommand(G4String("/process/em/verbose    0"));
+	UI->ApplyCommand(G4String("/process/had/verbose   0"));
+	UI->ApplyCommand(G4String("/process/eLoss/verbose 0"));
+
+	UI->ApplyCommand(G4String("/control/verbose  ") + (verbose > 0 ? "1" : "0"));
+	UI->ApplyCommand(G4String("/run/verbose      ") + (verbose > 1 ? "1" : "0"));
+	UI->ApplyCommand(G4String("/event/verbose    ") + (verbose > 2 ? "1" : "0"));
+	UI->ApplyCommand(G4String("/hits/verbose     ") + (verbose > 3 ? "1" : "0"));
+	UI->ApplyCommand(G4String("/tracking/verbose ") + (verbose > 4 ? "1" : "0"));
+	UI->ApplyCommand(G4String("/stepping/verbose ") + (verbose > 5 ? "1" : "0"));
+
+#ifdef G4MULTITHREADED
+	if (nthreads > 1)
+		runManager->SetNumberOfThreads(nthreads);
+#endif
+
+	runManager->SetUserInitialization(new DetectorConstruction(gdml_filename));
+	runManager->SetUserInitialization(factory.GetReferencePhysList(physics_list));
+	runManager->SetUserInitialization(new InitializationAction());
+
+	runManager->Initialize();
+
+	if (!runManager->ConfirmBeamOnCondition())
+		errx(1, "Geant4 is not fully initialized");
+
+	if (macro)
+		UI->ApplyCommand(G4String("/control/execute ") + macro);
+
+	if (nevents >= 0)
+		runManager->BeamOn(nevents);
+
+	if (interactive) {
+		enable_stdout();
+		for (;;) {
+			char command[1024];
+			fprintf(stdout, ">>> ");
+			if (!fgets(command, sizeof(command), stdin))
+				break;
+			fprintf(stderr, "command = '%s'\n", command);
+			UI->ApplyCommand(command);
+		}
+		if (!verbose)
+			disable_stdout();
+	}
+
+	delete runManager;
+	return 0;
+}

primary.cc

+#include "primary.h"
+
+#include <G4Event.hh>
+#include <G4SystemOfUnits.hh>
+#include <G4ParticleTable.hh>
+#include <G4PrimaryParticle.hh>
+#include <G4RandomDirection.hh>
+
+#include <err.h>
+
+static double primary_energy;
+static const char *primary_name;
+static const G4ParticleDefinition *primary;
+
+void set_primary_name(const char *name)
+{
+    primary_name = name;
+}
+
+void set_primary_energy(double E)
+{
+    primary_energy = E * GeV;
+}
+
+PrimaryGeneratorAction::PrimaryGeneratorAction()
+{
+    if (!(primary = G4ParticleTable::GetParticleTable()->FindParticle(primary_name)))
+        errx(1, "unknown particle type: %s", primary_name);
+}
+
+void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event)
+{
+	static const G4double time = 0.0;
+	static const G4ThreeVector position(0.0, 0.0, 0.0);
+
+	G4PrimaryVertex* vertex = new G4PrimaryVertex(position, time);
+	G4PrimaryParticle* particle = new G4PrimaryParticle(primary);
+
+	particle->SetKineticEnergy(primary_energy);
+	particle->SetMomentumDirection(G4RandomDirection());
+	particle->SetPolarization(0.0, 0.0, 0.0);
+	vertex->SetPrimary(particle);
+	event->AddPrimaryVertex(vertex);
+}

primary.h

+#ifndef _PRIMARY_H
+#define _PRIMARY_H
+
+#include <G4VUserPrimaryGeneratorAction.hh>
+
+void set_primary_name(const char*);
+void set_primary_energy(double);
+
+class G4Event;
+
+class PrimaryGeneratorAction final : public G4VUserPrimaryGeneratorAction {
+public:
+    PrimaryGeneratorAction();
+	void GeneratePrimaries(G4Event*) override;
+};
+
+#endif