diff --git a/.github/workflows/build-test.yaml b/.github/workflows/build-test.yaml index 1b4c768..be31e2c 100644 --- a/.github/workflows/build-test.yaml +++ b/.github/workflows/build-test.yaml @@ -15,10 +15,9 @@ on: pull_request: push: - branches: - [main] + branches: [main] - ## Paste this snippet into the workflow file to enable tmate debugging + ## Paste this snippet into the workflow file to enable tmate debugging # - name: Setup tmate session # uses: mxschmitt/action-tmate@v3 # if: @@ -208,6 +207,10 @@ jobs: command: | python -c "import geant4_python_application as g4; g4.install_datasets()" + - name: Run tests debug + run: | + python -m pytest -s -vv tests/test_analysis.py::test_sensitive + - name: Run tests run: | python -m pytest -vv --reruns 3 --reruns-delay 30 --only-rerun "(?i)http|timeout|connection|socket|resolve" diff --git a/pyproject.toml b/pyproject.toml index 765c332..cbf0618 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,6 +32,7 @@ classifiers = [ dependencies = [ "awkward", "vector", + "hepunits", "fsspec", "pyarrow", "awkward-pandas", diff --git a/src/geant4_application/include/geant4_application/Application.h b/src/geant4_application/include/geant4_application/Application.h index 17b8bb1..8fd3d46 100644 --- a/src/geant4_application/include/geant4_application/Application.h +++ b/src/geant4_application/include/geant4_application/Application.h @@ -49,7 +49,7 @@ class Application { void SetupAction(); void Initialize(); - py::list Run(const py::object& primaries); + std::vector Run(const py::object& primaries); bool IsSetup() const; bool IsInitialized() const; diff --git a/src/geant4_application/include/geant4_application/DataModel.h b/src/geant4_application/include/geant4_application/DataModel.h index f90755a..bb7ead0 100644 --- a/src/geant4_application/include/geant4_application/DataModel.h +++ b/src/geant4_application/include/geant4_application/DataModel.h @@ -178,7 +178,7 @@ void InsertEvent(const G4Event* event, Builders& builder); void InsertTrack(const G4Track* track, Builders& builder); void InsertStep(const G4Step* step, Builders& builder); -py::object SnapshotBuilder(Builders& builder); +py::object BuilderToObject(std::unique_ptr builders); namespace units { static constexpr auto energy = CLHEP::keV; diff --git a/src/geant4_application/include/geant4_application/EventAction.h b/src/geant4_application/include/geant4_application/EventAction.h index 22f637f..b819f0e 100644 --- a/src/geant4_application/include/geant4_application/EventAction.h +++ b/src/geant4_application/include/geant4_application/EventAction.h @@ -13,6 +13,8 @@ class EventAction : public G4UserEventAction { void BeginOfEventAction(const G4Event*) override; void EndOfEventAction(const G4Event*) override; + + static double sensitiveEnergy; }; }// namespace geant4_app diff --git a/src/geant4_application/include/geant4_application/RunAction.h b/src/geant4_application/include/geant4_application/RunAction.h index 9498740..3e1633f 100644 --- a/src/geant4_application/include/geant4_application/RunAction.h +++ b/src/geant4_application/include/geant4_application/RunAction.h @@ -1,6 +1,7 @@ #pragma once +#include #include #include @@ -20,12 +21,12 @@ class RunAction : public G4UserRunAction { /// Only one instance of RunAction is created for each thread. static data::Builders& GetBuilder(); - static std::unique_ptr GetContainer(); + static std::unique_ptr> GetContainer(); private: std::unique_ptr builder = nullptr; std::mutex mutex; - static std::unique_ptr container; + static std::unique_ptr> container; static std::vector> buildersToSnapshot; }; diff --git a/src/geant4_application/src/Application.cpp b/src/geant4_application/src/Application.cpp index beed4d8..e807553 100644 --- a/src/geant4_application/src/Application.cpp +++ b/src/geant4_application/src/Application.cpp @@ -32,10 +32,11 @@ Application::Application() { Application::~Application() = default; void Application::SetupRandomEngine() { - G4Random::setTheEngine(new CLHEP::RanecuEngine); + G4Random::setTheEngine(new CLHEP::MTwistEngine); if (randomSeed == 0) { randomSeed = std::random_device()(); } + cout << "seed set to: " << randomSeed << endl; G4Random::setTheSeed(randomSeed); } @@ -88,6 +89,10 @@ void Application::SetupManager(unsigned short nThreads) { delete G4VSteppingVerbose::GetInstance(); SteppingVerbose::SetInstance(new SteppingVerbose); + // https://geant4-forum.web.cern.ch/t/different-random-seeds-but-same-results/324/5 + // seed needs to be setup before the run manager is created + SetupRandomEngine(); + const auto runManagerType = nThreads > 0 ? G4RunManagerType::MTOnly : G4RunManagerType::SerialOnly; runManager = unique_ptr(G4RunManagerFactory::CreateRunManager(runManagerType)); if (nThreads > 0) { @@ -110,13 +115,11 @@ void Application::Initialize() { throw runtime_error("Application is already initialized"); } - SetupRandomEngine(); - runManager->Initialize(); isInitialized = true; } -py::list Application::Run(const py::object& primaries) { +vector Application::Run(const py::object& primaries) { if (!IsInitialized()) { Initialize(); } diff --git a/src/geant4_application/src/DataModel.cpp b/src/geant4_application/src/DataModel.cpp index f1c2dd9..b966cb8 100644 --- a/src/geant4_application/src/DataModel.cpp +++ b/src/geant4_application/src/DataModel.cpp @@ -429,82 +429,83 @@ py::object snapshot_builder(const T& builder) { return from_buffers(builder.form(), builder.length(), container); } -py::object SnapshotBuilder(Builders& builder) { +py::object BuilderToObject(unique_ptr builder) { + // this will automatically clear the builders after the pointer goes out of scope py::dict snapshot; - if (builder.fields.contains("run")) { - snapshot["run"] = snapshot_builder(builder.run); + if (builder->fields.contains("run")) { + snapshot["run"] = snapshot_builder(builder->run); } - if (builder.fields.contains("id")) { - snapshot["id"] = snapshot_builder(builder.id); + if (builder->fields.contains("id")) { + snapshot["id"] = snapshot_builder(builder->id); } - if (builder.fields.contains("primaries")) { - snapshot["primaries"] = snapshot_builder(builder.primaries); + if (builder->fields.contains("primaries")) { + snapshot["primaries"] = snapshot_builder(builder->primaries); } - if (builder.fields.contains("track_id")) { - snapshot["track_id"] = snapshot_builder(builder.track_id); + if (builder->fields.contains("track_id")) { + snapshot["track_id"] = snapshot_builder(builder->track_id); } - if (builder.fields.contains("track_parent_id")) { - snapshot["track_parent_id"] = snapshot_builder(builder.track_parent_id); + if (builder->fields.contains("track_parent_id")) { + snapshot["track_parent_id"] = snapshot_builder(builder->track_parent_id); } - if (builder.fields.contains("track_initial_energy")) { - snapshot["track_initial_energy"] = snapshot_builder(builder.track_initial_energy); + if (builder->fields.contains("track_initial_energy")) { + snapshot["track_initial_energy"] = snapshot_builder(builder->track_initial_energy); } - if (builder.fields.contains("track_initial_time")) { - snapshot["track_initial_time"] = snapshot_builder(builder.track_initial_time); + if (builder->fields.contains("track_initial_time")) { + snapshot["track_initial_time"] = snapshot_builder(builder->track_initial_time); } - if (builder.fields.contains("track_weight")) { - snapshot["track_weight"] = snapshot_builder(builder.track_weight); + if (builder->fields.contains("track_weight")) { + snapshot["track_weight"] = snapshot_builder(builder->track_weight); } - if (builder.fields.contains("track_initial_position")) { - snapshot["track_initial_position"] = snapshot_builder(builder.track_initial_position); + if (builder->fields.contains("track_initial_position")) { + snapshot["track_initial_position"] = snapshot_builder(builder->track_initial_position); } - if (builder.fields.contains("track_initial_momentum")) { - snapshot["track_initial_momentum"] = snapshot_builder(builder.track_initial_momentum); + if (builder->fields.contains("track_initial_momentum")) { + snapshot["track_initial_momentum"] = snapshot_builder(builder->track_initial_momentum); } - if (builder.fields.contains("track_particle")) { - snapshot["track_particle"] = snapshot_builder(builder.track_particle); + if (builder->fields.contains("track_particle")) { + snapshot["track_particle"] = snapshot_builder(builder->track_particle); } - if (builder.fields.contains("track_particle_type")) { - snapshot["track_particle_type"] = snapshot_builder(builder.track_particle_type); + if (builder->fields.contains("track_particle_type")) { + snapshot["track_particle_type"] = snapshot_builder(builder->track_particle_type); } - if (builder.fields.contains("track_creator_process")) { - snapshot["track_creator_process"] = snapshot_builder(builder.track_creator_process); + if (builder->fields.contains("track_creator_process")) { + snapshot["track_creator_process"] = snapshot_builder(builder->track_creator_process); } - if (builder.fields.contains("track_creator_process_type")) { - snapshot["track_creator_process_type"] = snapshot_builder(builder.track_creator_process_type); + if (builder->fields.contains("track_creator_process_type")) { + snapshot["track_creator_process_type"] = snapshot_builder(builder->track_creator_process_type); } - if (builder.fields.contains("track_children_ids")) { - snapshot["track_children_ids"] = snapshot_builder(builder.track_children_ids); + if (builder->fields.contains("track_children_ids")) { + snapshot["track_children_ids"] = snapshot_builder(builder->track_children_ids); } - if (builder.fields.contains("step_energy")) { - snapshot["step_energy"] = snapshot_builder(builder.step_energy); + if (builder->fields.contains("step_energy")) { + snapshot["step_energy"] = snapshot_builder(builder->step_energy); } - if (builder.fields.contains("step_time")) { - snapshot["step_time"] = snapshot_builder(builder.step_time); + if (builder->fields.contains("step_time")) { + snapshot["step_time"] = snapshot_builder(builder->step_time); } - if (builder.fields.contains("step_track_kinetic_energy")) { - snapshot["step_track_kinetic_energy"] = snapshot_builder(builder.step_track_kinetic_energy); + if (builder->fields.contains("step_track_kinetic_energy")) { + snapshot["step_track_kinetic_energy"] = snapshot_builder(builder->step_track_kinetic_energy); } - if (builder.fields.contains("step_process")) { - snapshot["step_process"] = snapshot_builder(builder.step_process); + if (builder->fields.contains("step_process")) { + snapshot["step_process"] = snapshot_builder(builder->step_process); } - if (builder.fields.contains("step_process_type")) { - snapshot["step_process_type"] = snapshot_builder(builder.step_process_type); + if (builder->fields.contains("step_process_type")) { + snapshot["step_process_type"] = snapshot_builder(builder->step_process_type); } - if (builder.fields.contains("step_volume")) { - snapshot["step_volume"] = snapshot_builder(builder.step_volume); + if (builder->fields.contains("step_volume")) { + snapshot["step_volume"] = snapshot_builder(builder->step_volume); } - if (builder.fields.contains("step_volume_post")) { - snapshot["step_volume_post"] = snapshot_builder(builder.step_volume_post); + if (builder->fields.contains("step_volume_post")) { + snapshot["step_volume_post"] = snapshot_builder(builder->step_volume_post); } - if (builder.fields.contains("step_nucleus")) { - snapshot["step_nucleus"] = snapshot_builder(builder.step_nucleus); + if (builder->fields.contains("step_nucleus")) { + snapshot["step_nucleus"] = snapshot_builder(builder->step_nucleus); } - if (builder.fields.contains("step_position")) { - snapshot["step_position"] = snapshot_builder(builder.step_position); + if (builder->fields.contains("step_position")) { + snapshot["step_position"] = snapshot_builder(builder->step_position); } - if (builder.fields.contains("step_momentum")) { - snapshot["step_momentum"] = snapshot_builder(builder.step_momentum); + if (builder->fields.contains("step_momentum")) { + snapshot["step_momentum"] = snapshot_builder(builder->step_momentum); } return snapshot; } diff --git a/src/geant4_application/src/EventAction.cpp b/src/geant4_application/src/EventAction.cpp index c821043..8090b22 100644 --- a/src/geant4_application/src/EventAction.cpp +++ b/src/geant4_application/src/EventAction.cpp @@ -11,10 +11,15 @@ using namespace geant4_app; EventAction::EventAction() : G4UserEventAction() {} void EventAction::BeginOfEventAction(const G4Event* event) { + sensitiveEnergy = 0; data::InsertEventBegin(event, RunAction::GetBuilder()); data::InsertEvent(event, RunAction::GetBuilder()); } void EventAction::EndOfEventAction(const G4Event* event) { data::InsertEventEnd(event, RunAction::GetBuilder()); + + cout << "END OF EVENT: event id: " << event->GetEventID() << " Sensitive energy: " << sensitiveEnergy << endl; } + +double EventAction::sensitiveEnergy = 0.0; diff --git a/src/geant4_application/src/PrimaryGeneratorAction.cpp b/src/geant4_application/src/PrimaryGeneratorAction.cpp index b365606..d1be609 100644 --- a/src/geant4_application/src/PrimaryGeneratorAction.cpp +++ b/src/geant4_application/src/PrimaryGeneratorAction.cpp @@ -16,34 +16,8 @@ using namespace geant4_app; PrimaryGeneratorAction::PrimaryGeneratorAction() : G4VUserPrimaryGeneratorAction() {} void PrimaryGeneratorAction::GeneratePrimaries(G4Event* event) { - if (!awkwardPrimaryEnergies.empty()) { - const double energy = awkwardPrimaryEnergies[event->GetEventID()]; - gun.SetParticleEnergy(energy * keV); - } - if (!awkwardPrimaryPositions.empty()) { - const auto& position = awkwardPrimaryPositions[event->GetEventID()]; - gun.SetParticlePosition(G4ThreeVector(position[0] * cm, position[1] * cm, position[2] * cm)); - } - if (!awkwardPrimaryDirections.empty()) { - const auto& direction = awkwardPrimaryDirections[event->GetEventID()]; - gun.SetParticleMomentumDirection(G4ThreeVector(direction[0], direction[1], direction[2])); - } - if (!awkwardPrimaryParticles.empty()) { - const auto& particleAwkward = awkwardPrimaryParticles[event->GetEventID()]; - auto* particle = G4ParticleTable::GetParticleTable()->FindParticle(particleAwkward); - if (particle == nullptr) { - throw runtime_error("PrimaryGeneratorAction::GeneratePrimaries - particle '" + particleAwkward + "' not found"); - } - gun.SetParticleDefinition(particle); - } - if (generatorType == "gun") { - gun.GeneratePrimaryVertex(event); - } else if (generatorType == "gps") { - gps.GeneratePrimaryVertex(event); - } else { - throw runtime_error("PrimaryGeneratorAction::GeneratePrimaries - generatorType must be 'gun', 'gps'"); - } + gun.GeneratePrimaryVertex(event); } void PrimaryGeneratorAction::SetGeneratorType(const string& type) { diff --git a/src/geant4_application/src/RunAction.cpp b/src/geant4_application/src/RunAction.cpp index b42cc57..13f2dcb 100644 --- a/src/geant4_application/src/RunAction.cpp +++ b/src/geant4_application/src/RunAction.cpp @@ -2,6 +2,8 @@ #include "geant4_application/RunAction.h" #include "geant4_application/SteppingVerbose.h" +#include + #include using namespace std; @@ -9,7 +11,7 @@ using namespace geant4_app; RunAction::RunAction() : G4UserRunAction() {} -void RunAction::BeginOfRunAction(const G4Run*) { +void RunAction::BeginOfRunAction(const G4Run* run) { lock_guard lock(mutex); builder = make_unique(Application::GetEventFields()); @@ -18,25 +20,33 @@ void RunAction::BeginOfRunAction(const G4Run*) { steppingVerbose->Initialize(); if (IsMaster()) { - container = make_unique(); + container = make_unique>(); + } + + cout << "RUN ID: " << run->GetRunID() << " RANDOM SEED: " << G4Random::getTheSeed() << endl; + + for (int i = 0; i < 10; ++i) { + cout << "RANDOM NUMBER: " << G4UniformRand() << endl; } } void RunAction::EndOfRunAction(const G4Run*) { lock_guard lock(mutex); + cout << "END OF RUN. Thread: " << G4Threading::G4GetThreadId() << endl; if (!isMaster || !G4Threading::IsMultithreadedApplication()) { buildersToSnapshot.push_back(std::move(builder)); + cout << "LENGTH: " << buildersToSnapshot.back()->run.length() << endl; } if (isMaster) { - for (auto& builderToSnapshot: buildersToSnapshot) { - auto data = SnapshotBuilder(*builderToSnapshot); - container->append(data); - builderToSnapshot = nullptr; + for (auto& builders: buildersToSnapshot) { + container->push_back(BuilderToObject(std::move(builders))); } buildersToSnapshot.clear(); } + + cout << "- END OF RUN. Thread: " << G4Threading::G4GetThreadId() << endl; } data::Builders& RunAction::GetBuilder() { @@ -44,9 +54,9 @@ data::Builders& RunAction::GetBuilder() { return *runAction->builder; } -unique_ptr RunAction::container = nullptr; +unique_ptr> RunAction::container = nullptr; -unique_ptr RunAction::GetContainer() { +unique_ptr> RunAction::GetContainer() { return std::move(RunAction::container); } diff --git a/src/geant4_application/src/SensitiveDetector.cpp b/src/geant4_application/src/SensitiveDetector.cpp index 5774e39..c8a2582 100644 --- a/src/geant4_application/src/SensitiveDetector.cpp +++ b/src/geant4_application/src/SensitiveDetector.cpp @@ -1,11 +1,15 @@ #include "geant4_application/SensitiveDetector.h" +#include "geant4_application/EventAction.h" using namespace std; using namespace geant4_app; SensitiveDetector::SensitiveDetector(const string& name) : G4VSensitiveDetector(name) {} -G4bool SensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) { return true; } +G4bool SensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) { + EventAction::sensitiveEnergy += step->GetTotalEnergyDeposit() / CLHEP::keV; + return true; +} void SensitiveDetector::Initialize(G4HCofThisEvent*) {} diff --git a/tests/test_analysis.py b/tests/test_analysis.py index 4e64f96..535158f 100644 --- a/tests/test_analysis.py +++ b/tests/test_analysis.py @@ -12,8 +12,7 @@ def test_events(): - with Application() as app: - app.seed = 1000 + with Application(seed=1000) as app: app.setup_detector(gdml=complexGdml) sensitive_volume = "gasVolume" @@ -42,31 +41,34 @@ def test_events(): def test_sensitive(): - with Application() as app: - app.seed = 1000 - app.setup_detector(gdml=complexGdml).setup_action() + with Application(seed=1234) as app: + app.setup_detector(gdml=complexGdml) - sensitive_volume = "gasVolume" - app.detector.sensitive_volumes = {sensitive_volume} + app.detector.sensitive_volumes = {"gasVolume"} + volume_logical = "gasVolume" app.initialize() - app.command("/gun/particle neutron") - app.command("/gun/energy 100 MeV") - app.command("/gun/direction 0 -1 0") - app.command("/gun/position 0 10 0 m") + app.command("/gun/particle gamma") + app.command("/gun/energy 10 keV") + app.command("/gun/direction 0 0 -1") + app.command("/gun/position 0 0 20 cm") - volumes = app.detector.physical_volumes_from_logical(sensitive_volume) + volumes = app.detector.physical_volumes_from_logical(volume_logical) assert len(volumes) == 1 volume = list(volumes)[0] print("volume: ", volume) - events = [] - total = 2 - while (n := np.sum([len(_events) for _events in events])) < total: - run_events = app.run(20) - run_events = run_events[run_events.energy_in_volume(volume) > 0] - events.append(run_events) - events = ak.concatenate(events) - # ak.to_parquet(events.hits(volume), "hits.parquet") + events = app.run(1000) + energy_in_sensitive = events.energy_in_volume(volume) + for i, energy in enumerate(energy_in_sensitive): + print(f"event {i} energy: {energy}") + + events = events[energy_in_sensitive > 0] + print("n: ", len(events)) + return + + assert app.seed == 1234 + assert len(events) == 5 hits = events.hits(volume) + np.isclose(hits.energy[0][0:5], [0.0401, 0.00271, 0.00391, 0.389, 0.204]) electrons = ak.concatenate([hit.electrons() for hit in hits]) diff --git a/tests/test_application.py b/tests/test_application.py index f43d1f0..377d0a9 100644 --- a/tests/test_application.py +++ b/tests/test_application.py @@ -47,15 +47,9 @@ def test_seed_single_thread(): # launch 100 events events = app.run(100) assert len(events) == 100 - reference_value = [ - 0.00000000e00, - 4.54930276e-01, - 2.25814551e-01, - 8.09506886e-03, - 1.22345221e00, - ] + reference_value = [0, 0.502, 0, 0.0952, 0.0199] energy = np.array(events.track.step.energy[0][0][0:5]) - assert np.allclose(energy, reference_value, atol=1e-5) + assert np.allclose(energy, reference_value, atol=1e-4) assert app.seed == 1100 diff --git a/tests/test_run.py b/tests/test_run.py index e58df82..630f510 100644 --- a/tests/test_run.py +++ b/tests/test_run.py @@ -4,6 +4,7 @@ import numpy as np import pytest import requests +from hepunits import MeV, keV from geant4_python_application import Application, basic_gdml @@ -34,3 +35,28 @@ def test_awkward_primaries(n_threads): event_primary_energy = ak.flatten(events.primaries.energy) assert np.allclose(event_primary_energy, primaries.energy) + + +@pytest.mark.parametrize("n_threads", [0]) +def test_reproduce(n_threads): + primaries = ak.Array( + [ + { + "particle": "neutron", + "energy": 100 * MeV / keV, + "direction": {"x": 0, "y": -1, "z": 0}, + "position": {"x": 0, "y": 10, "z": 0}, + } + for _ in range(500) + ] + ) + with Application(gdml=complex_gdml, n_threads=n_threads, seed=137) as app: + events = ak.Array([]) + batch_size = 100 + for i in range(0, len(primaries), batch_size): + events = ak.concatenate([events, app.run(primaries[i : i + batch_size])]) + + assert len(events) == len(primaries) + volume = app.detector.physical_volumes_from_logical("gasVolume") + events = events[events.energy_in_volume(volume) > 0] + assert len(events) == 7