From b75bdd799724ab5a9506f4aea39b03a7692717ce Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Thu, 9 Apr 2026 02:26:32 +0000 Subject: [PATCH 1/4] Add viztracks example for Opticks photon trajectory visualization Add examples/viztracks/ that converts GPU photon step records into G4VTrajectory objects for Geant4 interactive visualization. Uses KeepTheCurrentEvent to preserve the event, then injects trajectories from the Opticks record array in EndOfRunAction. --- CMakeLists.txt | 1 + examples/viztracks/CMakeLists.txt | 11 ++ examples/viztracks/OpticksTrajectory.h | 48 +++++ examples/viztracks/viztracks.cpp | 256 +++++++++++++++++++++++++ 4 files changed, 316 insertions(+) create mode 100644 examples/viztracks/CMakeLists.txt create mode 100644 examples/viztracks/OpticksTrajectory.h create mode 100644 examples/viztracks/viztracks.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index fe1b25c0d..0acce859e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -40,6 +40,7 @@ add_subdirectory(u4/tests) add_subdirectory(g4cx) add_subdirectory(g4cx/tests) add_subdirectory(src) +add_subdirectory(examples/viztracks) # Export configs include(CMakePackageConfigHelpers) diff --git a/examples/viztracks/CMakeLists.txt b/examples/viztracks/CMakeLists.txt new file mode 100644 index 000000000..4b141242c --- /dev/null +++ b/examples/viztracks/CMakeLists.txt @@ -0,0 +1,11 @@ +find_package(Geant4 REQUIRED ui_all vis_all) + +add_executable(viztracks viztracks.cpp) +target_link_libraries(viztracks gphox ${Geant4_LIBRARIES} G4CX) +target_include_directories(viztracks PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR} + $ + $ +) + +install(TARGETS viztracks RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/viztracks/OpticksTrajectory.h b/examples/viztracks/OpticksTrajectory.h new file mode 100644 index 000000000..36d38475b --- /dev/null +++ b/examples/viztracks/OpticksTrajectory.h @@ -0,0 +1,48 @@ +#pragma once + +#include + +#include "G4ThreeVector.hh" +#include "G4VTrajectory.hh" +#include "G4VTrajectoryPoint.hh" + +struct OpticksTrajectoryPoint : public G4VTrajectoryPoint +{ + G4ThreeVector fPosition; + + OpticksTrajectoryPoint(const G4ThreeVector &pos) : fPosition(pos) {} + ~OpticksTrajectoryPoint() override = default; + + const G4ThreeVector GetPosition() const override { return fPosition; } +}; + +struct OpticksTrajectory : public G4VTrajectory +{ + G4int fTrackID; + G4ThreeVector fInitialMomentum; + std::vector fPoints; + + OpticksTrajectory(G4int trackID, const G4ThreeVector &initialMomentum) + : fTrackID(trackID), fInitialMomentum(initialMomentum) + { + } + + ~OpticksTrajectory() override + { + for (auto *p : fPoints) + delete p; + } + + void AddPoint(const G4ThreeVector &pos) { fPoints.push_back(new OpticksTrajectoryPoint(pos)); } + + G4int GetTrackID() const override { return fTrackID; } + G4int GetParentID() const override { return 0; } + G4String GetParticleName() const override { return "opticalphoton"; } + G4double GetCharge() const override { return 0.; } + G4int GetPDGEncoding() const override { return 0; } + G4ThreeVector GetInitialMomentum() const override { return fInitialMomentum; } + G4int GetPointEntries() const override { return G4int(fPoints.size()); } + G4VTrajectoryPoint *GetPoint(G4int i) const override { return fPoints[i]; } + void AppendStep(const G4Step *) override {} + void MergeTrajectory(G4VTrajectory *) override {} +}; diff --git a/examples/viztracks/viztracks.cpp b/examples/viztracks/viztracks.cpp new file mode 100644 index 000000000..e49f574f6 --- /dev/null +++ b/examples/viztracks/viztracks.cpp @@ -0,0 +1,256 @@ +/** + * viztracks — Visualize Opticks GPU photon trajectories in Geant4 + * + * Runs the same Geant4 + Opticks GPU simulation as GPURaytrace, then + * converts the photon step records into G4VTrajectory objects and + * injects them into the G4Event for interactive visualization. + * + * Requires: OPTICKS_EVENT_MODE=DebugHeavy (or DebugLite) and + * OPTICKS_MAX_SLOT=M1 to gather record arrays from GPU. + * + * Usage: + * OPTICKS_EVENT_MODE=DebugHeavy OPTICKS_MAX_SLOT=M1 \ + * viztracks -g geometry.gdml -m run.mac -s 42 -i + */ + +#include + +#include + +#include "FTFP_BERT.hh" +#include "G4OpticalPhysics.hh" +#include "G4VModularPhysicsList.hh" + +#include "G4UIExecutive.hh" +#include "G4UImanager.hh" +#include "G4VisExecutive.hh" + +#include "sysrap/OPTICKS_LOG.hh" + +#include "src/GPURaytrace.h" + +#include "G4EventManager.hh" +#include "G4Run.hh" +#include "G4RunManager.hh" +#include "G4RunManagerFactory.hh" +#include "G4TrajectoryContainer.hh" +#include "G4VUserActionInitialization.hh" + +#include "sysrap/NPFold.h" +#include "sysrap/SEventConfig.hh" +#include "sysrap/sphoton.h" + +#include "OpticksTrajectory.h" + +using namespace std; + +struct VizRunAction : G4UserRunAction +{ + EventAction *fEventAction; + + VizRunAction(EventAction *eventAction) : fEventAction(eventAction) {} + + void BeginOfRunAction(const G4Run *) override {} + + void EndOfRunAction(const G4Run *run) override + { + if (!G4Threading::IsMasterThread()) + return; + + G4CXOpticks *gx = G4CXOpticks::Get(); + + auto start = std::chrono::high_resolution_clock::now(); + gx->simulate(0, false); + cudaDeviceSynchronize(); + auto end = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed = end - start; + cout << "Simulation time: " << elapsed.count() << " seconds" << endl; + + SEvt *sev = SEvt::Get_EGPU(); + unsigned int num_hits = sev->GetNumHit(0); + + cout << "Opticks: NumGensteps: " << sev->GetNumGenstepFromGenstep(0) << endl; + cout << "Opticks: NumPhotons: " << sev->GetNumPhotonCollected(0) << endl; + cout << "Opticks: NumHits: " << num_hits << endl; + cout << "Geant4: NumHits: " << fEventAction->GetTotalG4Hits() << endl; + + // Build trajectories from the photon step record array + const NP *record = sev->topfold->get("record"); + if (!record || record->shape.size() < 2) + { + cerr << "No record array — set OPTICKS_EVENT_MODE=DebugHeavy and OPTICKS_MAX_SLOT=M1" << endl; + return; + } + + int num_photons = record->shape[0]; + int max_steps = record->shape[1]; + const sphoton *steps = (const sphoton *)record->bytes(); + + auto *trajContainer = new G4TrajectoryContainer(); + + for (int j = 0; j < num_photons; j++) + { + const sphoton &first = steps[j * max_steps]; + auto *traj = new OpticksTrajectory( + j, G4ThreeVector(first.mom.x, first.mom.y, first.mom.z)); + + for (int i = 0; i < max_steps; i++) + { + const sphoton &s = steps[j * max_steps + i]; + if (s.flagmask == 0) + break; + traj->AddPoint(G4ThreeVector(s.pos.x, s.pos.y, s.pos.z)); + } + + if (traj->GetPointEntries() >= 2) + trajContainer->push_back(traj); + else + delete traj; + } + + // Inject trajectories into the last kept event + auto *eventVector = run->GetEventVector(); + G4Event *event = + (eventVector && !eventVector->empty()) + ? const_cast(eventVector->back()) + : nullptr; + + if (event) + { + event->SetTrajectoryContainer(trajContainer); + cout << "Injected " << trajContainer->size() + << " Opticks photon trajectories into G4Event" << endl; + } + else + { + delete trajContainer; + cerr << "No kept G4Event available for trajectory injection" << endl; + } + } +}; + +struct VizEventAction : public EventAction +{ + using EventAction::EventAction; + + void EndOfEventAction(const G4Event *event) override + { + EventAction::EndOfEventAction(event); + // Keep the event alive so EndOfRunAction can inject trajectories + G4EventManager::GetEventManager()->KeepTheCurrentEvent(); + } +}; + +struct VizActionInitialization : public G4VUserActionInitialization +{ + G4App *fG4App; + VizEventAction *fVizEventAction; + VizRunAction *fVizRunAction; + + VizActionInitialization(G4App *app) + : G4VUserActionInitialization(), fG4App(app), + fVizEventAction(new VizEventAction(app->sev)), + fVizRunAction(new VizRunAction(fVizEventAction)) + { + } + + void BuildForMaster() const override + { + SetUserAction(fVizRunAction); + } + + void Build() const override + { + SetUserAction(fG4App->prim_gen_); + SetUserAction(fVizRunAction); + SetUserAction(fVizEventAction); + SetUserAction(fG4App->tracking_); + SetUserAction(fG4App->stepping_); + } +}; + +int main(int argc, char **argv) +{ + OPTICKS_LOG(argc, argv); + + argparse::ArgumentParser program("viztracks", "0.0.0"); + + string gdml_file, macro_name; + bool interactive; + + program.add_argument("-g", "--gdml") + .help("path to GDML file") + .default_value(string("geom.gdml")) + .nargs(1) + .store_into(gdml_file); + + program.add_argument("-m", "--macro") + .help("path to G4 macro") + .default_value(string("run.mac")) + .nargs(1) + .store_into(macro_name); + + program.add_argument("-i", "--interactive") + .help("open an interactive viewer (required to see trajectories)") + .flag() + .store_into(interactive); + + program.add_argument("-s", "--seed").help("fixed random seed (default: time-based)").scan<'i', long>(); + + try + { + program.parse_args(argc, argv); + } + catch (const exception &err) + { + cerr << err.what() << endl; + cerr << program; + exit(EXIT_FAILURE); + } + + long seed; + if (program.is_used("--seed")) + { + seed = program.get("--seed"); + } + else + { + seed = static_cast(time(nullptr)); + } + CLHEP::HepRandom::setTheSeed(seed); + G4cout << "Random seed set to: " << seed << G4endl; + + G4VModularPhysicsList *physics = new FTFP_BERT; + physics->RegisterPhysics(new G4OpticalPhysics); + + auto *run_mgr = G4RunManagerFactory::CreateRunManager(); + run_mgr->SetUserInitialization(physics); + + G4App *g4app = new G4App(gdml_file); + + VizActionInitialization *actionInit = new VizActionInitialization(g4app); + run_mgr->SetUserInitialization(actionInit); + run_mgr->SetUserInitialization(g4app->det_cons_); + + G4UIExecutive *uix = nullptr; + G4VisManager *vis = nullptr; + + if (interactive) + { + uix = new G4UIExecutive(argc, argv); + vis = new G4VisExecutive; + vis->Initialize(); + } + + G4UImanager *ui = G4UImanager::GetUIpointer(); + ui->ApplyCommand("/control/execute " + macro_name); + + if (interactive) + { + uix->SessionStart(); + } + + delete uix; + + return EXIT_SUCCESS; +} From 79f68eaf385ee47d5da8c9f8c530ef53ba6307c3 Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Thu, 9 Apr 2026 13:37:55 +0000 Subject: [PATCH 2/4] Add AlignGpuCpuTracks example for photon-by-photon GPU vs G4 comparison Aligned validation using precooked curand Philox sequences so G4 consumes the same random numbers as the GPU. Includes: - StandAloneGeant4Validation with ShimG4Op* for RNG alignment - generate_precooked_rng.cu for Philox sequence generation - compare_aligned.py for flag/position/chi2 comparison - raindrop configs for testing Tested on raindrop geometry: 10000/10000 photons match (100%). --- CMakeLists.txt | 1 + examples/AlignGpuCpuTracks/CMakeLists.txt | 15 + .../StandAloneGeant4Validation.cpp | 167 +++++ .../StandAloneGeant4Validation.h | 578 ++++++++++++++++++ examples/AlignGpuCpuTracks/compare_aligned.py | 207 +++++++ .../generate_precooked_rng.cu | 113 ++++ examples/AlignGpuCpuTracks/raindrop.json | 30 + .../AlignGpuCpuTracks/raindrop_sphere.json | 30 + 8 files changed, 1141 insertions(+) create mode 100644 examples/AlignGpuCpuTracks/CMakeLists.txt create mode 100644 examples/AlignGpuCpuTracks/StandAloneGeant4Validation.cpp create mode 100644 examples/AlignGpuCpuTracks/StandAloneGeant4Validation.h create mode 100644 examples/AlignGpuCpuTracks/compare_aligned.py create mode 100644 examples/AlignGpuCpuTracks/generate_precooked_rng.cu create mode 100644 examples/AlignGpuCpuTracks/raindrop.json create mode 100644 examples/AlignGpuCpuTracks/raindrop_sphere.json diff --git a/CMakeLists.txt b/CMakeLists.txt index 0acce859e..5b2a50193 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -41,6 +41,7 @@ add_subdirectory(g4cx) add_subdirectory(g4cx/tests) add_subdirectory(src) add_subdirectory(examples/viztracks) +add_subdirectory(examples/AlignGpuCpuTracks) # Export configs include(CMakePackageConfigHelpers) diff --git a/examples/AlignGpuCpuTracks/CMakeLists.txt b/examples/AlignGpuCpuTracks/CMakeLists.txt new file mode 100644 index 000000000..c9bf9cfe7 --- /dev/null +++ b/examples/AlignGpuCpuTracks/CMakeLists.txt @@ -0,0 +1,15 @@ +find_package(Geant4 REQUIRED ui_all vis_all) + +add_executable(StandAloneGeant4Validation + StandAloneGeant4Validation.cpp + StandAloneGeant4Validation.h +) +target_link_libraries(StandAloneGeant4Validation gphox ${Geant4_LIBRARIES} U4) +target_compile_definitions(StandAloneGeant4Validation PRIVATE WITH_INSTRUMENTED_DEBUG) +target_include_directories(StandAloneGeant4Validation PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR} + $ + $ +) + +install(TARGETS StandAloneGeant4Validation RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.cpp b/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.cpp new file mode 100644 index 000000000..259c0d6d8 --- /dev/null +++ b/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.cpp @@ -0,0 +1,167 @@ +#include +#include + +#include + +#include "FTFP_BERT.hh" +#include "G4OpticalPhysics.hh" +#include "G4MTRunManager.hh" +#include "G4RunManager.hh" +#include "G4VModularPhysicsList.hh" +#include "G4UImanager.hh" + +#include "G4OpticalParameters.hh" + +#include "StandAloneGeant4Validation.h" +#include "src/config.h" + +using namespace std; + +int main(int argc, char **argv) +{ + argparse::ArgumentParser program("StandAloneGeant4Validation", "0.0.0"); + + string gdml_file, config_name; + int num_threads = 0; + + program.add_argument("-g", "--gdml") + .help("path to GDML file") + .default_value(string("geom.gdml")) + .nargs(1) + .store_into(gdml_file); + + program.add_argument("-c", "--config") + .help("the name of a config file") + .default_value(string("dev")) + .nargs(1) + .store_into(config_name); + + program.add_argument("-s", "--seed") + .help("fixed random seed (default: time-based)") + .scan<'i', long>(); + + program.add_argument("-t", "--threads") + .help("number of threads (0=sequential, default: hardware concurrency)") + .default_value(-1) + .scan<'i', int>() + .store_into(num_threads); + + program.add_argument("--aligned") + .help("enable photon-by-photon aligned comparison with GPU (forces sequential)") + .default_value(false) + .implicit_value(true); + + try + { + program.parse_args(argc, argv); + } + catch (const exception &err) + { + cerr << err.what() << endl; + cerr << program; + exit(EXIT_FAILURE); + } + + long seed; + if (program.is_used("--seed")) + seed = program.get("--seed"); + else + seed = static_cast(time(nullptr)); + + bool aligned = program.get("--aligned"); + + gphox::Config cfg(config_name); + int total_photons = cfg.torch.numphoton; + + // Aligned mode forces sequential (U4Random is single-threaded) + if (aligned) + num_threads = 0; + + // Determine threading mode + bool use_mt = (num_threads != 0); + if (num_threads < 0) + num_threads = std::thread::hardware_concurrency(); + if (num_threads < 1) + num_threads = 1; + + // In MT mode: split photons across events, one event per thread-batch + // In sequential mode: one event with all photons (original behavior) + int num_events, photons_per_event; + if (use_mt) + { + num_events = num_threads * 4; // 4 events per thread for load balancing + photons_per_event = (total_photons + num_events - 1) / num_events; + // Adjust num_events so we don't overshoot + num_events = (total_photons + photons_per_event - 1) / photons_per_event; + } + else + { + num_events = 1; + photons_per_event = total_photons; + } + + int actual_photons = num_events * photons_per_event; + + G4cout << "Random seed set to: " << seed << G4endl; + G4cout << "G4: " << total_photons << " photons, " + << num_events << " events x " << photons_per_event << " photons/event" + << " (" << actual_photons << " actual)" + << (use_mt ? ", " + to_string(num_threads) + " threads" : ", sequential") + << G4endl; + + HitAccumulator accumulator; + PhotonFateAccumulator fate; + + if (aligned) + fate.Resize(total_photons); + + G4VModularPhysicsList *physics = new FTFP_BERT; + if (aligned) + physics->RegisterPhysics(new AlignedOpticalPhysics); + else + physics->RegisterPhysics(new G4OpticalPhysics); + + // Use exponential WLS time profile (default is delta = zero delay) + G4OpticalParameters::Instance()->SetWLSTimeProfile("exponential"); + + if (use_mt) + { + auto *run_mgr = new G4MTRunManager; + run_mgr->SetNumberOfThreads(num_threads); + run_mgr->SetUserInitialization(physics); + run_mgr->SetUserInitialization(new G4OnlyDetectorConstruction(gdml_file, &accumulator)); + run_mgr->SetUserInitialization( + new G4OnlyActionInitialization(cfg, &accumulator, &fate, photons_per_event, num_events, aligned)); + run_mgr->Initialize(); + + CLHEP::HepRandom::setTheSeed(seed); + + G4cout << "G4: Starting MT run with " << num_events << " events..." << G4endl; + run_mgr->BeamOn(num_events); + + delete run_mgr; + } + else + { + G4RunManager run_mgr; + run_mgr.SetUserInitialization(physics); + run_mgr.SetUserInitialization(new G4OnlyDetectorConstruction(gdml_file, &accumulator)); + run_mgr.SetUserInitialization( + new G4OnlyActionInitialization(cfg, &accumulator, &fate, photons_per_event, num_events, aligned)); + + if (aligned) + { + G4cout << "G4: Aligned mode — creating U4Random" << G4endl; + U4Random::Create(); + } + + run_mgr.Initialize(); + + CLHEP::HepRandom::setTheSeed(seed); + + G4cout << "G4: Starting sequential run..." << G4endl; + run_mgr.BeamOn(num_events); + } + + return EXIT_SUCCESS; +} diff --git a/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.h b/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.h new file mode 100644 index 000000000..540fea9f0 --- /dev/null +++ b/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.h @@ -0,0 +1,578 @@ +#pragma once + +#include +#include +#include +#include + +#include "G4Event.hh" +#include "G4GDMLParser.hh" +#include "G4THitsCollection.hh" +#include "G4VHit.hh" +#include "G4OpticalPhoton.hh" +#include "G4PhysicalConstants.hh" +#include "G4PrimaryParticle.hh" +#include "G4PrimaryVertex.hh" +#include "G4Run.hh" +#include "G4SDManager.hh" +#include "G4SystemOfUnits.hh" +#include "G4ThreeVector.hh" +#include "G4Track.hh" +#include "G4TrackStatus.hh" +#include "G4UserEventAction.hh" +#include "G4UserRunAction.hh" +#include "G4UserSteppingAction.hh" +#include "G4UserTrackingAction.hh" +#include "G4VPhysicalVolume.hh" +#include "G4VUserActionInitialization.hh" +#include "G4VUserDetectorConstruction.hh" +#include "G4VUserPrimaryGeneratorAction.hh" +#include "G4OpBoundaryProcess.hh" +#include "G4ProcessManager.hh" +#include "G4VPhysicsConstructor.hh" +#include "G4OpWLS.hh" + +#include "ShimG4OpAbsorption.hh" +#include "ShimG4OpRayleigh.hh" +#include "U4Random.hh" + +#include "sysrap/NP.hh" +#include "sysrap/sphoton.h" + +#include "src/config.h" +#include "src/torch.h" + +// ---- Global hit accumulator (thread-safe) ---- + +struct HitAccumulator +{ + std::mutex mtx; + std::vector hits; + + void AddHits(const std::vector &event_hits) + { + std::lock_guard lock(mtx); + hits.insert(hits.end(), event_hits.begin(), event_hits.end()); + } + + void Save(const char *filename) + { + std::lock_guard lock(mtx); + G4int num_hits = hits.size(); + NP *arr = NP::Make(num_hits, 4, 4); + for (int i = 0; i < num_hits; i++) + { + float *data = reinterpret_cast(&hits[i]); + std::copy(data, data + 16, arr->values() + i * 16); + } + arr->save(filename); + delete arr; + G4cout << "G4: Saved " << num_hits << " total hits to " << filename << G4endl; + } +}; + +// ---- Sensitive Detector: collects optical photon hits per event ---- + +struct G4PhotonHit : public G4VHit +{ + G4PhotonHit() = default; + + G4PhotonHit(G4double energy, G4double time, G4ThreeVector position, + G4ThreeVector direction, G4ThreeVector polarization) + : photon() + { + photon.pos = {static_cast(position.x()), + static_cast(position.y()), + static_cast(position.z())}; + photon.time = time; + photon.mom = {static_cast(direction.x()), + static_cast(direction.y()), + static_cast(direction.z())}; + photon.pol = {static_cast(polarization.x()), + static_cast(polarization.y()), + static_cast(polarization.z())}; + photon.wavelength = h_Planck * c_light / (energy * CLHEP::eV); + } + + void Print() override { G4cout << photon << G4endl; } + + sphoton photon; +}; + +using G4PhotonHitsCollection = G4THitsCollection; + +struct G4PhotonSD : public G4VSensitiveDetector +{ + HitAccumulator *accumulator; + + G4PhotonSD(G4String name, HitAccumulator *acc) + : G4VSensitiveDetector(name), accumulator(acc) + { + G4String HCname = name + "_HC"; + collectionName.insert(HCname); + } + + void Initialize(G4HCofThisEvent *hce) override + { + fHitsCollection = new G4PhotonHitsCollection(SensitiveDetectorName, collectionName[0]); + if (fHCID < 0) + fHCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); + hce->AddHitsCollection(fHCID, fHitsCollection); + } + + G4bool ProcessHits(G4Step *aStep, G4TouchableHistory *) override + { + G4Track *track = aStep->GetTrack(); + if (track->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) + return false; + + G4PhotonHit *hit = new G4PhotonHit( + track->GetTotalEnergy(), + track->GetGlobalTime(), + aStep->GetPostStepPoint()->GetPosition(), + aStep->GetPostStepPoint()->GetMomentumDirection(), + aStep->GetPostStepPoint()->GetPolarization()); + + fHitsCollection->insert(hit); + track->SetTrackStatus(fStopAndKill); + return true; + } + + void EndOfEvent(G4HCofThisEvent *) override + { + G4int num_hits = fHitsCollection->entries(); + + std::vector event_hits; + event_hits.reserve(num_hits); + for (G4PhotonHit *hit : *fHitsCollection->GetVector()) + event_hits.push_back(hit->photon); + + accumulator->AddHits(event_hits); + } + + private: + G4PhotonHitsCollection *fHitsCollection = nullptr; + G4int fHCID = -1; +}; + +// ---- Detector Construction: loads GDML, attaches SD ---- + +struct G4OnlyDetectorConstruction : G4VUserDetectorConstruction +{ + G4OnlyDetectorConstruction(std::filesystem::path gdml_file, HitAccumulator *acc) + : gdml_file_(gdml_file), accumulator_(acc) {} + + G4VPhysicalVolume *Construct() override + { + parser_.Read(gdml_file_.string(), false); + return parser_.GetWorldVolume(); + } + + void ConstructSDandField() override + { + G4SDManager *SDman = G4SDManager::GetSDMpointer(); + const G4GDMLAuxMapType *auxmap = parser_.GetAuxMap(); + + for (auto const &[logVol, listType] : *auxmap) + { + for (auto const &auxtype : listType) + { + if (auxtype.type == "SensDet") + { + G4String name = logVol->GetName() + "_" + auxtype.value; + G4cout << "G4: Attaching SD to " << logVol->GetName() << G4endl; + G4PhotonSD *sd = new G4PhotonSD(name, accumulator_); + SDman->AddNewDetector(sd); + logVol->SetSensitiveDetector(sd); + } + } + } + } + + private: + std::filesystem::path gdml_file_; + G4GDMLParser parser_; + HitAccumulator *accumulator_; +}; + +// ---- Primary Generator: distributes photons across events ---- + +struct G4OnlyPrimaryGenerator : G4VUserPrimaryGeneratorAction +{ + gphox::Config cfg; + int photons_per_event; + + G4OnlyPrimaryGenerator(const gphox::Config &cfg, int photons_per_event) + : cfg(cfg), photons_per_event(photons_per_event) {} + + void GeneratePrimaries(G4Event *event) override + { + int eventID = event->GetEventID(); + + // Generate photons for this event's batch using event-specific seed offset + storch t = cfg.torch; + t.numphoton = photons_per_event; + std::vector sphotons = generate_photons(t, photons_per_event, eventID); + + for (const sphoton &p : sphotons) + { + G4ThreeVector position(p.pos.x, p.pos.y, p.pos.z); + G4ThreeVector direction(p.mom.x, p.mom.y, p.mom.z); + G4ThreeVector polarization(p.pol.x, p.pol.y, p.pol.z); + G4double wavelength_nm = p.wavelength; + + G4PrimaryVertex *vertex = new G4PrimaryVertex(position, p.time); + G4double energy = h_Planck * c_light / (wavelength_nm * nm); + + G4PrimaryParticle *particle = new G4PrimaryParticle(G4OpticalPhoton::Definition()); + particle->SetKineticEnergy(energy); + particle->SetMomentumDirection(direction); + particle->SetPolarization(polarization); + + vertex->SetPrimary(particle); + event->AddPrimaryVertex(vertex); + } + } +}; + +// ---- Photon fate accumulator: tracks ALL photon final states ---- + +struct PhotonFateAccumulator +{ + std::mutex mtx; + std::vector photons; + bool indexed = false; // true for aligned mode: store by photon index + + // Opticks flag enum values + static constexpr unsigned TORCH = 0x0004; + static constexpr unsigned BULK_ABSORB = 0x0008; + static constexpr unsigned BULK_REEMIT = 0x0010; + static constexpr unsigned BULK_SCATTER = 0x0020; + static constexpr unsigned SURFACE_DETECT = 0x0040; + static constexpr unsigned SURFACE_ABSORB = 0x0080; + static constexpr unsigned SURFACE_DREFLECT = 0x0100; + static constexpr unsigned SURFACE_SREFLECT = 0x0200; + static constexpr unsigned BOUNDARY_REFLECT = 0x0400; + static constexpr unsigned BOUNDARY_TRANSMIT= 0x0800; + static constexpr unsigned MISS = 0x8000; + + void Resize(int n) + { + photons.resize(n); + indexed = true; + } + + void Set(int idx, const sphoton& p) + { + if (idx >= 0 && idx < (int)photons.size()) + photons[idx] = p; + } + + void Add(const sphoton& p) + { + std::lock_guard lock(mtx); + photons.push_back(p); + } + + void Save(const char* filename) + { + std::lock_guard lock(mtx); + int n = photons.size(); + NP* arr = NP::Make(n, 4, 4); + for (int i = 0; i < n; i++) + { + float* data = reinterpret_cast(&photons[i]); + std::copy(data, data + 16, arr->values() + i * 16); + } + arr->save(filename); + delete arr; + G4cout << "G4: Saved " << n << " photon fates to " << filename << G4endl; + } +}; + +// ---- Stepping Action: tracks photon fates with opticks-compatible flags ---- + +struct G4OnlySteppingAction : G4UserSteppingAction +{ + PhotonFateAccumulator* fate; + bool aligned; + std::map proc_death_counts; + std::map boundary_status_counts; + std::mutex count_mtx; + + G4OnlySteppingAction(PhotonFateAccumulator* f, bool aligned_ = false) + : fate(f), aligned(aligned_) {} + + ~G4OnlySteppingAction() + { + std::lock_guard lock(count_mtx); + if (!proc_death_counts.empty()) + { + G4cout << "\nG4: Photon death process summary:" << G4endl; + for (auto& [name, count] : proc_death_counts) + G4cout << " " << name << ": " << count << G4endl; + } + if (!boundary_status_counts.empty()) + { + G4cout << "\nG4: OpBoundary status counts (all steps):" << G4endl; + const char* bnames[] = { + "Undefined","Transmission","FresnelRefraction","FresnelReflection", + "TotalInternalReflection","LambertianReflection","LobeReflection", + "SpikeReflection","BackScattering","Absorption","Detection", + "NotAtBoundary","SameMaterial","StepTooSmall","NoRINDEX", + "PolishedLumirrorAirReflection","PolishedLumirrorGlueReflection", + "PolishedAirReflection","PolishedTeflonAirReflection", + "PolishedTiOAirReflection","PolishedTyvekAirReflection", + "PolishedVM2000AirReflection","PolishedVM2000GlueReflection", + "EtchedLumirrorAirReflection","EtchedLumirrorGlueReflection", + "EtchedAirReflection","EtchedTeflonAirReflection", + "EtchedTiOAirReflection","EtchedTyvekAirReflection", + "EtchedVM2000AirReflection","EtchedVM2000GlueReflection", + "GroundLumirrorAirReflection","GroundLumirrorGlueReflection", + "GroundAirReflection","GroundTeflonAirReflection", + "GroundTiOAirReflection","GroundTyvekAirReflection", + "GroundVM2000AirReflection","GroundVM2000GlueReflection", + "Dichroic","CoatedDielectricReflection","CoatedDielectricRefraction", + "CoatedDielectricFrustratedTransmission" + }; + for (auto& [st, count] : boundary_status_counts) + { + const char* nm = (st >= 0 && st < 43) ? bnames[st] : "?"; + G4cout << " " << nm << "(" << st << "): " << count << G4endl; + } + } + } + + void UserSteppingAction(const G4Step* aStep) override + { + G4Track* track = aStep->GetTrack(); + if (track->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) + return; + + G4StepPoint* post = aStep->GetPostStepPoint(); + G4TrackStatus status = track->GetTrackStatus(); + + // Find the OpBoundary process to get its status (for ALL steps) + G4OpBoundaryProcess* boundary = nullptr; + G4ProcessManager* pm = track->GetDefinition()->GetProcessManager(); + for (int i = 0; i < pm->GetPostStepProcessVector()->entries(); i++) + { + G4VProcess* p = (*pm->GetPostStepProcessVector())[i]; + boundary = dynamic_cast(p); + if (boundary) break; + } + + G4OpBoundaryProcessStatus bStatus = boundary ? boundary->GetStatus() : Undefined; + + // Count boundary status for ALL steps + if (boundary && bStatus != NotAtBoundary && bStatus != Undefined && bStatus != StepTooSmall) + { + std::lock_guard lock(count_mtx); + boundary_status_counts[int(bStatus)]++; + } + + // Only record photon state when the photon is about to die + if (status != fStopAndKill && status != fStopButAlive) + return; + + // Identify the process + const G4VProcess* proc = post->GetProcessDefinedStep(); + G4String procName = proc ? proc->GetProcessName() : "Unknown"; + + // Build detailed key for counting + std::string key = procName; + if (procName == "OpBoundary" && boundary) + key += "(" + std::to_string(int(bStatus)) + ")"; + key += (status == fStopAndKill ? "/Kill" : "/Alive"); + + { + std::lock_guard lock(count_mtx); + proc_death_counts[key]++; + } + + // Map to opticks flag + unsigned flag = 0; + + if (procName == "OpAbsorption") + { + flag = PhotonFateAccumulator::BULK_ABSORB; + } + else if (procName == "OpWLS") + { + flag = PhotonFateAccumulator::BULK_REEMIT; + } + else if (procName == "OpBoundary" && boundary) + { + switch (bStatus) + { + case Detection: flag = PhotonFateAccumulator::SURFACE_DETECT; break; + case Absorption: flag = PhotonFateAccumulator::SURFACE_ABSORB; break; + case FresnelReflection: + case TotalInternalReflection: + flag = PhotonFateAccumulator::BOUNDARY_REFLECT; break; + case FresnelRefraction: flag = PhotonFateAccumulator::BOUNDARY_TRANSMIT; break; + case LambertianReflection: + case LobeReflection: flag = PhotonFateAccumulator::SURFACE_DREFLECT; break; + case SpikeReflection: flag = PhotonFateAccumulator::SURFACE_SREFLECT; break; + case BackScattering: flag = PhotonFateAccumulator::SURFACE_DREFLECT; break; + default: flag = PhotonFateAccumulator::SURFACE_ABSORB; break; + } + } + else if (procName == "Transportation") + { + // Check if an SD killed this photon (SURFACE_DETECT) + G4StepPoint* pre = aStep->GetPreStepPoint(); + G4VPhysicalVolume* preVol = pre->GetPhysicalVolume(); + G4VPhysicalVolume* postVol = post->GetPhysicalVolume(); + G4LogicalVolume* preLog = preVol ? preVol->GetLogicalVolume() : nullptr; + G4LogicalVolume* postLog = postVol ? postVol->GetLogicalVolume() : nullptr; + bool sd_pre = preLog && preLog->GetSensitiveDetector(); + bool sd_post = postLog && postLog->GetSensitiveDetector(); + if (sd_pre || sd_post) + flag = PhotonFateAccumulator::SURFACE_DETECT; + else + flag = PhotonFateAccumulator::BOUNDARY_TRANSMIT; + } + + if (flag == 0) flag = PhotonFateAccumulator::MISS; // catch-all + + // Build sphoton with the final state + G4ThreeVector pos = post->GetPosition(); + G4ThreeVector mom = post->GetMomentumDirection(); + G4ThreeVector pol = post->GetPolarization(); + G4double time = post->GetGlobalTime(); + G4double energy = post->GetTotalEnergy(); + + sphoton p = {}; + p.pos = { float(pos.x()), float(pos.y()), float(pos.z()) }; + p.time = float(time); + p.mom = { float(mom.x()), float(mom.y()), float(mom.z()) }; + p.pol = { float(pol.x()), float(pol.y()), float(pol.z()) }; + p.wavelength = (energy > 0) ? float(h_Planck * c_light / (energy * CLHEP::eV)) : 0.f; + + p.orient_boundary_flag = flag & 0xFFFF; + p.flagmask = flag; + + if (aligned && fate->indexed) + { + int photon_idx = track->GetTrackID() - 1; // G4 trackIDs are 1-based + fate->Set(photon_idx, p); + } + else + { + fate->Add(p); + } + } +}; + +// ---- Tracking Action: per-photon RNG sync for aligned mode ---- + +struct G4OnlyTrackingAction : G4UserTrackingAction +{ + void PreUserTrackingAction(const G4Track* track) override + { + if (track->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) + return; + int photon_idx = track->GetTrackID() - 1; // G4 trackIDs are 1-based + U4Random::SetSequenceIndex(photon_idx); + } + + void PostUserTrackingAction(const G4Track* track) override + { + if (track->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) + return; + U4Random::SetSequenceIndex(-1); + } +}; + +// ---- AlignedOpticalPhysics: uses Shim processes for precise RNILL matching ---- + +struct AlignedOpticalPhysics : G4VPhysicsConstructor +{ + AlignedOpticalPhysics() : G4VPhysicsConstructor("AlignedOptical") {} + void ConstructParticle() override {} + void ConstructProcess() override + { + auto* pm = G4OpticalPhoton::OpticalPhoton()->GetProcessManager(); + pm->AddDiscreteProcess(new ShimG4OpAbsorption()); + pm->AddDiscreteProcess(new ShimG4OpRayleigh()); + pm->AddDiscreteProcess(new G4OpBoundaryProcess()); + pm->AddDiscreteProcess(new G4OpWLS()); + } +}; + +// ---- Event Action: reports per-event progress ---- + +struct G4OnlyEventAction : G4UserEventAction +{ + int total_events; + + G4OnlyEventAction(int total_events) : total_events(total_events) {} + + void EndOfEventAction(const G4Event *event) override + { + int id = event->GetEventID(); + if (id == 0 || (id + 1) % 10 == 0 || id + 1 == total_events) + G4cout << "G4: Event " << id + 1 << "/" << total_events << G4endl; + } +}; + +// ---- Run Action: saves merged hits at end ---- + +struct G4OnlyRunAction : G4UserRunAction +{ + HitAccumulator *accumulator; + PhotonFateAccumulator *fate; + + G4OnlyRunAction(HitAccumulator *acc, PhotonFateAccumulator *f = nullptr) + : accumulator(acc), fate(f) {} + + void EndOfRunAction(const G4Run *) override + { + if (G4Threading::IsMasterThread() || !G4Threading::IsMultithreadedApplication()) + { + G4cout << "G4: Total accumulated hits: " << accumulator->hits.size() << G4endl; + accumulator->Save("g4_hits.npy"); + if (fate) + { + G4cout << "G4: Total photon fates: " << fate->photons.size() << G4endl; + fate->Save("g4_photon.npy"); + } + } + } +}; + +// ---- Action Initialization (required for MT) ---- + +struct G4OnlyActionInitialization : G4VUserActionInitialization +{ + gphox::Config cfg; + HitAccumulator *accumulator; + PhotonFateAccumulator *fate; + int photons_per_event; + int num_events; + bool aligned; + + G4OnlyActionInitialization(const gphox::Config &cfg, HitAccumulator *acc, + PhotonFateAccumulator *f, + int photons_per_event, int num_events, + bool aligned_ = false) + : cfg(cfg), accumulator(acc), fate(f), + photons_per_event(photons_per_event), + num_events(num_events), aligned(aligned_) {} + + void BuildForMaster() const override + { + SetUserAction(new G4OnlyRunAction(accumulator, fate)); + } + + void Build() const override + { + SetUserAction(new G4OnlyPrimaryGenerator(cfg, photons_per_event)); + SetUserAction(new G4OnlyEventAction(num_events)); + SetUserAction(new G4OnlyRunAction(accumulator, fate)); + SetUserAction(new G4OnlySteppingAction(fate, aligned)); + if (aligned) + SetUserAction(new G4OnlyTrackingAction()); + } +}; diff --git a/examples/AlignGpuCpuTracks/compare_aligned.py b/examples/AlignGpuCpuTracks/compare_aligned.py new file mode 100644 index 000000000..b044ce6f2 --- /dev/null +++ b/examples/AlignGpuCpuTracks/compare_aligned.py @@ -0,0 +1,207 @@ +#!/usr/bin/env python3 +""" +compare_aligned.py - Photon-by-photon comparison of GPU vs G4 aligned simulations. + +Usage: + python compare_aligned.py + +Performs: + 1. Per-photon flag comparison (exact match rate) + 2. Position comparison at multiple thresholds + 3. Chi-squared test on flag distributions (gold-standard validation metric) + 4. Glancing-angle photon identification (normal sign ambiguity) + 5. Divergent photon listing +""" +import sys +import numpy as np + +FLAG_NAMES = { + 0x0004: "TORCH", 0x0008: "BULK_ABSORB", 0x0010: "BULK_REEMIT", + 0x0020: "BULK_SCATTER", 0x0040: "SURFACE_DETECT", 0x0080: "SURFACE_ABSORB", + 0x0100: "SURFACE_DREFLECT", 0x0200: "SURFACE_SREFLECT", + 0x0400: "BOUNDARY_REFLECT", 0x0800: "BOUNDARY_TRANSMIT", 0x8000: "MISS", +} + +def flag_name(f): + return FLAG_NAMES.get(f, f"0x{f:04x}") + +def extract_flag(photon): + """Extract flag from q3.x (orient_boundary_flag) - lower 16 bits.""" + q3 = photon.view(np.uint32).reshape(-1, 4, 4) + return q3[:, 3, 0] & 0xFFFF + +def chi2_flag_distribution(gpu_flags, g4_flags): + """ + Chi-squared comparison of flag distributions. + + Compares the frequency of each flag value between GPU and G4. + This is the opticks gold-standard validation metric. + + Returns (chi2, ndof, flags_used, gpu_counts, g4_counts). + """ + all_flags = sorted(set(gpu_flags) | set(g4_flags)) + gpu_counts = np.array([(gpu_flags == f).sum() for f in all_flags], dtype=float) + g4_counts = np.array([(g4_flags == f).sum() for f in all_flags], dtype=float) + + total = gpu_counts + g4_counts + mask = total > 0 + gpu_c = gpu_counts[mask] + g4_c = g4_counts[mask] + tot = total[mask] + flags_used = [f for f, m in zip(all_flags, mask) if m] + + n_gpu = gpu_c.sum() + n_g4 = g4_c.sum() + expected_gpu = tot * n_gpu / (n_gpu + n_g4) + expected_g4 = tot * n_g4 / (n_gpu + n_g4) + + chi2 = 0.0 + for i in range(len(flags_used)): + if expected_gpu[i] > 0: + chi2 += (gpu_c[i] - expected_gpu[i])**2 / expected_gpu[i] + if expected_g4[i] > 0: + chi2 += (g4_c[i] - expected_g4[i])**2 / expected_g4[i] + + ndof = max(len(flags_used) - 1, 1) + return chi2, ndof, flags_used, gpu_c, g4_c + +def identify_glancing(gpu, g4): + """ + Identify glancing-angle photons where the normal sign ambiguity + causes momentum negation between GPU and G4. + + At glancing incidence cos(theta) ~ 0, float32 vs float64 can produce + opposite normal signs, reflecting the photon in the opposite direction. + These photons have matching flags but very different positions. + + Returns boolean mask of glancing photons. + """ + gpu_mom = gpu[:, 1, :3] + g4_mom = g4[:, 1, :3] + + # Normalize momenta (should already be unit vectors, but be safe) + gpu_norm = np.linalg.norm(gpu_mom, axis=1, keepdims=True) + g4_norm = np.linalg.norm(g4_mom, axis=1, keepdims=True) + gpu_norm[gpu_norm == 0] = 1 + g4_norm[g4_norm == 0] = 1 + + gpu_hat = gpu_mom / gpu_norm + g4_hat = g4_mom / g4_norm + + # Dot product of momentum directions: -1 = fully negated (normal flip) + mom_dot = np.sum(gpu_hat * g4_hat, axis=1) + + # Glancing: momentum vectors are nearly anti-parallel (dot ~ -1) + glancing = mom_dot < -0.5 + return glancing, mom_dot + +def main(): + if len(sys.argv) < 3: + print(f"Usage: {sys.argv[0]} ") + sys.exit(1) + + gpu = np.load(sys.argv[1]) + g4 = np.load(sys.argv[2]) + + print(f"GPU shape: {gpu.shape}") + print(f"G4 shape: {g4.shape}") + + n = min(len(gpu), len(g4)) + gpu = gpu[:n] + g4 = g4[:n] + + gpu_flags = extract_flag(gpu) + g4_flags = extract_flag(g4) + + # ---- 1. Per-photon flag comparison ---- + match = gpu_flags == g4_flags + n_match = match.sum() + n_diff = n - n_match + print(f"\n{'='*60}") + print(f"FLAG COMPARISON ({n} photons)") + print(f"{'='*60}") + print(f" Matching: {n_match} ({100*n_match/n:.1f}%)") + print(f" Differ: {n_diff} ({100*n_diff/n:.1f}%)") + + # ---- 2. Position comparison ---- + gpu_pos = gpu[:, 0, :3] + g4_pos = g4[:, 0, :3] + pos_diff = np.linalg.norm(gpu_pos - g4_pos, axis=1) + zero_g4 = np.all(g4_pos == 0, axis=1) + + valid = ~zero_g4 + n_valid = valid.sum() + print(f"\n{'='*60}") + print(f"POSITION COMPARISON ({n_valid} valid, {zero_g4.sum()} unrecorded)") + print(f"{'='*60}") + if n_valid > 0: + vdiff = pos_diff[valid] + print(f" Mean dist: {vdiff.mean():.4f} mm") + print(f" Max dist: {vdiff.max():.4f} mm") + print(f" < 0.01 mm: {(vdiff < 0.01).sum()} ({100*(vdiff < 0.01).sum()/n_valid:.1f}%)") + print(f" < 0.1 mm: {(vdiff < 0.1).sum()} ({100*(vdiff < 0.1).sum()/n_valid:.1f}%)") + print(f" < 1.0 mm: {(vdiff < 1.0).sum()} ({100*(vdiff < 1.0).sum()/n_valid:.1f}%)") + + # ---- 3. Chi-squared test on flag distributions ---- + print(f"\n{'='*60}") + print(f"CHI-SQUARED TEST (flag distribution)") + print(f"{'='*60}") + + chi2_val, ndof, flags_used, gpu_c, g4_c = chi2_flag_distribution(gpu_flags, g4_flags) + + print(f" {'Flag':<20s} {'GPU':>8s} {'G4':>8s} {'Diff':>8s}") + print(f" {'-'*20} {'-'*8} {'-'*8} {'-'*8}") + for i, f in enumerate(flags_used): + diff = int(gpu_c[i] - g4_c[i]) + sign = "+" if diff > 0 else "" + print(f" {flag_name(f):<20s} {int(gpu_c[i]):>8d} {int(g4_c[i]):>8d} {sign}{diff:>7d}") + + deviant_frac = 100 * n_diff / n if n > 0 else 0 + print(f"\n chi2/ndof = {chi2_val:.2f}/{ndof} = {chi2_val/ndof:.2f}") + print(f" deviant fraction: {deviant_frac:.2f}% ({n_diff}/{n})") + + # ---- 4. Glancing-angle analysis ---- + print(f"\n{'='*60}") + print(f"GLANCING-ANGLE ANALYSIS (normal sign ambiguity)") + print(f"{'='*60}") + + glancing, mom_dot = identify_glancing(gpu, g4) + n_glancing = glancing.sum() + + # Among matching-flag photons, how many are glancing with large pos diff? + match_glancing = match & glancing + match_large_pos = match & (pos_diff > 1.0) + match_glancing_large = match & glancing & (pos_diff > 1.0) + + print(f" Glancing photons (mom dot < -0.5): {n_glancing}") + print(f" Matching flag + pos diff > 1mm: {match_large_pos.sum()}") + print(f" Of those, glancing: {match_glancing_large.sum()}") + if match_large_pos.sum() > 0: + frac = 100 * match_glancing_large.sum() / match_large_pos.sum() + print(f" Fraction explained by glancing: {frac:.0f}%") + + # Position stats excluding glancing photons + non_glancing_match = match & ~glancing & valid + if non_glancing_match.sum() > 0: + ng_diff = pos_diff[non_glancing_match] + print(f"\n Position (matching, non-glancing, {non_glancing_match.sum()} photons):") + print(f" Max dist: {ng_diff.max():.6f} mm") + print(f" Mean dist: {ng_diff.mean():.6f} mm") + print(f" < 0.01 mm: {(ng_diff < 0.01).sum()} ({100*(ng_diff < 0.01).sum()/non_glancing_match.sum():.1f}%)") + + # ---- 5. Divergent photon listing ---- + if n_diff > 0: + div_idx = np.where(~match)[0] + print(f"\n{'='*60}") + print(f"DIVERGENT PHOTONS (first 10 of {n_diff})") + print(f"{'='*60}") + for i in div_idx[:10]: + gf = flag_name(gpu_flags[i]) + cf = flag_name(g4_flags[i]) + gp = gpu_pos[i] + cp = g4_pos[i] + print(f" [{i:5d}] GPU: {gf:20s} pos=({gp[0]:8.2f},{gp[1]:8.2f},{gp[2]:8.2f})") + print(f" G4: {cf:20s} pos=({cp[0]:8.2f},{cp[1]:8.2f},{cp[2]:8.2f})") + +if __name__ == "__main__": + main() diff --git a/examples/AlignGpuCpuTracks/generate_precooked_rng.cu b/examples/AlignGpuCpuTracks/generate_precooked_rng.cu new file mode 100644 index 000000000..79cb2910d --- /dev/null +++ b/examples/AlignGpuCpuTracks/generate_precooked_rng.cu @@ -0,0 +1,113 @@ +/** +generate_precooked_rng.cu +========================== + +Generates precooked curand Philox sequences for U4Random aligned mode. +Each photon gets its own random stream matching the GPU simulation. + +Build: + nvcc -o generate_precooked_rng tools/generate_precooked_rng.cu \ + -I. -I/opt/eic-opticks/include/eic-opticks -lcurand -std=c++17 + +Usage: + ./generate_precooked_rng [num_photons] [num_randoms_per_photon] + Defaults: 100000 photons, 256 randoms each (nj=16, nk=16) + +Output: + ~/.opticks/precooked/QSimTest/rng_sequence/ + rng_sequence_f_ni_nj_nk_tranche/ + rng_sequence_f_ni_nj_nk_ioffset000000.npy +*/ + +#include +#include +#include +#include +#include + +#include +#include "sysrap/NP.hh" + +__global__ void generate_sequences(float* out, unsigned ni, unsigned nv, unsigned id_offset) +{ + unsigned idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx >= ni) return; + + unsigned photon_idx = id_offset + idx; + + // Match GPU simulation: curand_init(seed=0, subsequence=photon_idx, offset=0) + curandStatePhilox4_32_10_t rng; + curand_init(0ULL, (unsigned long long)photon_idx, 0ULL, &rng); + + float* row = out + idx * nv; + for (unsigned j = 0; j < nv; j++) + row[j] = curand_uniform(&rng); +} + +static void mkdirp(const char* path) +{ + char tmp[1024]; + snprintf(tmp, sizeof(tmp), "%s", path); + for (char* p = tmp + 1; *p; p++) + { + if (*p == '/') { *p = 0; mkdir(tmp, 0755); *p = '/'; } + } + mkdir(tmp, 0755); +} + +int main(int argc, char** argv) +{ + unsigned ni = 100000; + unsigned nj = 16; + unsigned nk = 16; + + if (argc > 1) ni = atoi(argv[1]); + if (argc > 2) + { + unsigned total = atoi(argv[2]); + nj = 1; nk = total; + for (unsigned f = 2; f * f <= total; f++) + { + if (total % f == 0 && f <= 64) { nj = f; nk = total / f; } + } + } + + unsigned nv = nj * nk; + printf("Generating precooked curand Philox sequences:\n"); + printf(" photons: %u, randoms/photon: %u (nj=%u, nk=%u), memory: %.1f MB\n", + ni, nv, nj, nk, (double)ni * nv * sizeof(float) / (1024 * 1024)); + + const char* home = getenv("HOME"); + char dirpath[512], filename[256], fullpath[768]; + + snprintf(dirpath, sizeof(dirpath), + "%s/.opticks/precooked/QSimTest/rng_sequence/rng_sequence_f_ni%u_nj%u_nk%u_tranche%u", + home, ni, nj, nk, ni); + mkdirp(dirpath); + + snprintf(filename, sizeof(filename), + "rng_sequence_f_ni%u_nj%u_nk%u_ioffset%06u.npy", ni, nj, nk, 0); + snprintf(fullpath, sizeof(fullpath), "%s/%s", dirpath, filename); + + float* d_out = nullptr; + cudaMalloc(&d_out, (size_t)ni * nv * sizeof(float)); + + unsigned threads = 256; + unsigned blocks = (ni + threads - 1) / threads; + generate_sequences<<>>(d_out, ni, nv, 0); + cudaDeviceSynchronize(); + + cudaError_t err = cudaGetLastError(); + if (err != cudaSuccess) { fprintf(stderr, "CUDA error: %s\n", cudaGetErrorString(err)); return 1; } + + NP* seq = NP::Make(ni, nj, nk); + cudaMemcpy(seq->values(), d_out, (size_t)ni * nv * sizeof(float), cudaMemcpyDeviceToHost); + cudaFree(d_out); + + seq->save(fullpath); + printf("Saved: %s\n", fullpath); + printf("Set OPTICKS_RANDOM_SEQPATH=%s\n", fullpath); + + delete seq; + return 0; +} diff --git a/examples/AlignGpuCpuTracks/raindrop.json b/examples/AlignGpuCpuTracks/raindrop.json new file mode 100644 index 000000000..f8cc035f9 --- /dev/null +++ b/examples/AlignGpuCpuTracks/raindrop.json @@ -0,0 +1,30 @@ +{ + "torch": { + "gentype": "TORCH", + "trackid": 0, + "matline": 0, + "numphoton": 10000, + + "pos": [0.0, 0.0, 0.0], + "time": 0.0, + + "mom": [0.0, 0.0, 1.0], + "weight": 0.0, + + "pol": [1.0, 0.0, 0.0], + "wavelength": 420.0, + + "zenith": [0.0, 1.0], + "azimuth": [0.0, 1.0], + + "radius": 15.0, + "distance": 0.0, + "mode": 255, + "type": "disc" + }, + + "event": { + "mode": "DebugHeavy", + "maxslot": 1000000 + } +} diff --git a/examples/AlignGpuCpuTracks/raindrop_sphere.json b/examples/AlignGpuCpuTracks/raindrop_sphere.json new file mode 100644 index 000000000..31071a978 --- /dev/null +++ b/examples/AlignGpuCpuTracks/raindrop_sphere.json @@ -0,0 +1,30 @@ +{ + "torch": { + "gentype": "TORCH", + "trackid": 0, + "matline": 0, + "numphoton": 10000, + + "pos": [0.0, 0.0, 0.0], + "time": 0.0, + + "mom": [0.0, 0.0, 1.0], + "weight": 0.0, + + "pol": [1.0, 0.0, 0.0], + "wavelength": 420.0, + + "zenith": [0.0, 1.0], + "azimuth": [0.0, 1.0], + + "radius": 0.0, + "distance": 0.0, + "mode": 255, + "type": "sphere" + }, + + "event": { + "mode": "DebugHeavy", + "maxslot": 1000000 + } +} From c7edb2ad6d4c52063f3e46087bf3ee9e7a1bda66 Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Fri, 10 Apr 2026 01:08:18 +0000 Subject: [PATCH 3/4] Add G4 step recording and VTK track visualization - StandAloneGeant4Validation: add StepRecordAccumulator that saves per-step photon positions as g4_record.npy (N, 32, 4, 4), matching the Opticks record.npy format for photon-by-photon track comparison - Add plot_divergent_tracks.py: VTK offscreen renderer that overlays GPU (red) and G4 (blue) photon tracks on GDML wireframe geometry - Add apex.json config for apex.gdml torch photon simulation --- .../StandAloneGeant4Validation.cpp | 9 +- .../StandAloneGeant4Validation.h | 103 +++++++- examples/AlignGpuCpuTracks/apex.json | 30 +++ .../plot_divergent_tracks.py | 224 ++++++++++++++++++ 4 files changed, 355 insertions(+), 11 deletions(-) create mode 100644 examples/AlignGpuCpuTracks/apex.json create mode 100644 examples/AlignGpuCpuTracks/plot_divergent_tracks.py diff --git a/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.cpp b/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.cpp index 259c0d6d8..a4589bcbe 100644 --- a/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.cpp +++ b/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.cpp @@ -111,9 +111,13 @@ int main(int argc, char **argv) HitAccumulator accumulator; PhotonFateAccumulator fate; + StepRecordAccumulator *record = nullptr; if (aligned) + { fate.Resize(total_photons); + record = new StepRecordAccumulator(total_photons, 32); + } G4VModularPhysicsList *physics = new FTFP_BERT; if (aligned) @@ -131,7 +135,7 @@ int main(int argc, char **argv) run_mgr->SetUserInitialization(physics); run_mgr->SetUserInitialization(new G4OnlyDetectorConstruction(gdml_file, &accumulator)); run_mgr->SetUserInitialization( - new G4OnlyActionInitialization(cfg, &accumulator, &fate, photons_per_event, num_events, aligned)); + new G4OnlyActionInitialization(cfg, &accumulator, &fate, photons_per_event, num_events, aligned, record)); run_mgr->Initialize(); CLHEP::HepRandom::setTheSeed(seed); @@ -147,7 +151,7 @@ int main(int argc, char **argv) run_mgr.SetUserInitialization(physics); run_mgr.SetUserInitialization(new G4OnlyDetectorConstruction(gdml_file, &accumulator)); run_mgr.SetUserInitialization( - new G4OnlyActionInitialization(cfg, &accumulator, &fate, photons_per_event, num_events, aligned)); + new G4OnlyActionInitialization(cfg, &accumulator, &fate, photons_per_event, num_events, aligned, record)); if (aligned) { @@ -163,5 +167,6 @@ int main(int argc, char **argv) run_mgr.BeamOn(num_events); } + delete record; return EXIT_SUCCESS; } diff --git a/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.h b/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.h index 540fea9f0..9f4d6a957 100644 --- a/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.h +++ b/examples/AlignGpuCpuTracks/StandAloneGeant4Validation.h @@ -290,18 +290,69 @@ struct PhotonFateAccumulator } }; +// ---- Step Record: saves per-step photon state as (num_photons, max_steps, 4, 4) ---- + +struct StepRecordAccumulator +{ + int max_photons; + int max_steps; + std::vector records; // flat array: [photon_idx * max_steps + step_idx] + std::vector step_counts; // current step count per photon + + StepRecordAccumulator(int np, int ms) : max_photons(np), max_steps(ms), + records(np * ms), step_counts(np, 0) + { + memset(records.data(), 0, records.size() * sizeof(sphoton)); + } + + void RecordStep(int photon_idx, const G4StepPoint* point, unsigned flag) + { + if (photon_idx < 0 || photon_idx >= max_photons) return; + int step = step_counts[photon_idx]; + if (step >= max_steps) return; + + G4ThreeVector pos = point->GetPosition(); + G4ThreeVector mom = point->GetMomentumDirection(); + G4ThreeVector pol = point->GetPolarization(); + G4double energy = point->GetTotalEnergy(); + + sphoton& s = records[photon_idx * max_steps + step]; + s.pos = { float(pos.x()), float(pos.y()), float(pos.z()) }; + s.time = float(point->GetGlobalTime()); + s.mom = { float(mom.x()), float(mom.y()), float(mom.z()) }; + s.pol = { float(pol.x()), float(pol.y()), float(pol.z()) }; + s.wavelength = (energy > 0) ? float(h_Planck * c_light / (energy * CLHEP::eV)) : 0.f; + s.orient_boundary_flag = flag & 0xFFFF; + s.flagmask = flag; + + step_counts[photon_idx]++; + } + + void Save(const char* filename) + { + NP* arr = NP::Make(max_photons, max_steps, 4, 4); + memcpy(arr->values(), records.data(), records.size() * sizeof(sphoton)); + arr->save(filename); + delete arr; + G4cout << "G4: Saved step records (" << max_photons << " x " << max_steps + << ") to " << filename << G4endl; + } +}; + // ---- Stepping Action: tracks photon fates with opticks-compatible flags ---- struct G4OnlySteppingAction : G4UserSteppingAction { PhotonFateAccumulator* fate; + StepRecordAccumulator* record; bool aligned; std::map proc_death_counts; std::map boundary_status_counts; std::mutex count_mtx; - G4OnlySteppingAction(PhotonFateAccumulator* f, bool aligned_ = false) - : fate(f), aligned(aligned_) {} + G4OnlySteppingAction(PhotonFateAccumulator* f, bool aligned_ = false, + StepRecordAccumulator* rec = nullptr) + : fate(f), record(rec), aligned(aligned_) {} ~G4OnlySteppingAction() { @@ -349,9 +400,21 @@ struct G4OnlySteppingAction : G4UserSteppingAction if (track->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) return; + if (track->GetCurrentStepNumber() > 1000) + track->SetTrackStatus(fStopAndKill); + G4StepPoint* post = aStep->GetPostStepPoint(); G4TrackStatus status = track->GetTrackStatus(); + // Record every step for aligned photon-by-photon comparison + if (record) + { + int photon_idx = track->GetTrackID() - 1; + if (track->GetCurrentStepNumber() == 1) + record->RecordStep(photon_idx, aStep->GetPreStepPoint(), 0x0004); // TORCH + record->RecordStep(photon_idx, post, 0x0001); // placeholder flag + } + // Find the OpBoundary process to get its status (for ALL steps) G4OpBoundaryProcess* boundary = nullptr; G4ProcessManager* pm = track->GetDefinition()->GetProcessManager(); @@ -436,6 +499,22 @@ struct G4OnlySteppingAction : G4UserSteppingAction if (flag == 0) flag = PhotonFateAccumulator::MISS; // catch-all + // Update the last recorded step with the death flag + if (record) + { + int photon_idx = track->GetTrackID() - 1; + if (photon_idx >= 0 && photon_idx < record->max_photons) + { + int step = record->step_counts[photon_idx]; + if (step > 0) + { + sphoton& last = record->records[photon_idx * record->max_steps + step - 1]; + last.orient_boundary_flag = flag & 0xFFFF; + last.flagmask = flag; + } + } + } + // Build sphoton with the final state G4ThreeVector pos = post->GetPosition(); G4ThreeVector mom = post->GetMomentumDirection(); @@ -523,9 +602,11 @@ struct G4OnlyRunAction : G4UserRunAction { HitAccumulator *accumulator; PhotonFateAccumulator *fate; + StepRecordAccumulator *record; - G4OnlyRunAction(HitAccumulator *acc, PhotonFateAccumulator *f = nullptr) - : accumulator(acc), fate(f) {} + G4OnlyRunAction(HitAccumulator *acc, PhotonFateAccumulator *f = nullptr, + StepRecordAccumulator *rec = nullptr) + : accumulator(acc), fate(f), record(rec) {} void EndOfRunAction(const G4Run *) override { @@ -538,6 +619,8 @@ struct G4OnlyRunAction : G4UserRunAction G4cout << "G4: Total photon fates: " << fate->photons.size() << G4endl; fate->Save("g4_photon.npy"); } + if (record) + record->Save("g4_record.npy"); } } }; @@ -549,6 +632,7 @@ struct G4OnlyActionInitialization : G4VUserActionInitialization gphox::Config cfg; HitAccumulator *accumulator; PhotonFateAccumulator *fate; + StepRecordAccumulator *record; int photons_per_event; int num_events; bool aligned; @@ -556,22 +640,23 @@ struct G4OnlyActionInitialization : G4VUserActionInitialization G4OnlyActionInitialization(const gphox::Config &cfg, HitAccumulator *acc, PhotonFateAccumulator *f, int photons_per_event, int num_events, - bool aligned_ = false) - : cfg(cfg), accumulator(acc), fate(f), + bool aligned_ = false, + StepRecordAccumulator *rec = nullptr) + : cfg(cfg), accumulator(acc), fate(f), record(rec), photons_per_event(photons_per_event), num_events(num_events), aligned(aligned_) {} void BuildForMaster() const override { - SetUserAction(new G4OnlyRunAction(accumulator, fate)); + SetUserAction(new G4OnlyRunAction(accumulator, fate, record)); } void Build() const override { SetUserAction(new G4OnlyPrimaryGenerator(cfg, photons_per_event)); SetUserAction(new G4OnlyEventAction(num_events)); - SetUserAction(new G4OnlyRunAction(accumulator, fate)); - SetUserAction(new G4OnlySteppingAction(fate, aligned)); + SetUserAction(new G4OnlyRunAction(accumulator, fate, record)); + SetUserAction(new G4OnlySteppingAction(fate, aligned, record)); if (aligned) SetUserAction(new G4OnlyTrackingAction()); } diff --git a/examples/AlignGpuCpuTracks/apex.json b/examples/AlignGpuCpuTracks/apex.json new file mode 100644 index 000000000..31071a978 --- /dev/null +++ b/examples/AlignGpuCpuTracks/apex.json @@ -0,0 +1,30 @@ +{ + "torch": { + "gentype": "TORCH", + "trackid": 0, + "matline": 0, + "numphoton": 10000, + + "pos": [0.0, 0.0, 0.0], + "time": 0.0, + + "mom": [0.0, 0.0, 1.0], + "weight": 0.0, + + "pol": [1.0, 0.0, 0.0], + "wavelength": 420.0, + + "zenith": [0.0, 1.0], + "azimuth": [0.0, 1.0], + + "radius": 0.0, + "distance": 0.0, + "mode": 255, + "type": "sphere" + }, + + "event": { + "mode": "DebugHeavy", + "maxslot": 1000000 + } +} diff --git a/examples/AlignGpuCpuTracks/plot_divergent_tracks.py b/examples/AlignGpuCpuTracks/plot_divergent_tracks.py new file mode 100644 index 000000000..a44448033 --- /dev/null +++ b/examples/AlignGpuCpuTracks/plot_divergent_tracks.py @@ -0,0 +1,224 @@ +#!/usr/bin/env python3 +""" +plot_divergent_tracks.py - VTK offscreen render of divergent GPU vs G4 photon tracks + +Usage: + python plot_divergent_tracks.py [num_tracks] [output.png] + +Renders the detector geometry (wireframe) with divergent photon tracks overlaid: + Red = GPU tracks + Blue = G4 tracks +""" +import sys +import numpy as np +import vtk +import xml.etree.ElementTree as ET + + +def parse_gdml_boxes(gdml_path): + """Parse GDML for box solids and physvol placements.""" + tree = ET.parse(gdml_path) + root = tree.getroot() + + solids = {} + for box in root.iter('box'): + solids[box.get('name')] = (float(box.get('x')), float(box.get('y')), float(box.get('z'))) + + physvols = {} + for pv in root.iter('physvol'): + name = pv.get('name', '') + volref = pv.find('volumeref') + posref = pv.find('position') + if volref is not None and posref is not None: + physvols[name] = (volref.get('ref'), + float(posref.get('x', '0')), + float(posref.get('y', '0')), + float(posref.get('z', '0'))) + + volumes = {} + for vol in root.iter('volume'): + sref = vol.find('solidref') + if sref is not None: + volumes[vol.get('name')] = sref.get('ref') + + return solids, physvols, volumes + + +COLORS_MAP = { + 'pTPlayer': (0.6, 0, 0.8, 0.6), + 'pTPsubstrate': (0.2, 0.2, 0.9, 0.25), + 'BlueWLS': (0, 0.7, 0.9, 0.25), + 'ReflectiveFoilBack': (1, 0.5, 0, 0.6), + 'ReflectiveFoilEdge': (0.9, 0.5, 0, 0.4), + 'SiPMs': (0, 0.7, 0, 0.5), +} + + +def get_volume_color(name): + for key, col in COLORS_MAP.items(): + if key in name: + return col + return (0.5, 0.5, 0.5, 0.1) + + +def add_geometry(renderer, solids, physvols, volumes): + """Add GDML box volumes as wireframes.""" + for pvname, (vname, px, py, pz) in physvols.items(): + sname = volumes.get(vname) + if sname and sname in solids: + sx, sy, sz = solids[sname] + cube = vtk.vtkCubeSource() + cube.SetXLength(sx); cube.SetYLength(sy); cube.SetZLength(sz) + cube.SetCenter(px, py, pz); cube.Update() + mapper = vtk.vtkPolyDataMapper() + mapper.SetInputData(cube.GetOutput()) + actor = vtk.vtkActor() + actor.SetMapper(mapper) + r, g, b, a = get_volume_color(pvname) + actor.GetProperty().SetColor(r, g, b) + actor.GetProperty().SetOpacity(a) + actor.GetProperty().SetRepresentationToWireframe() + actor.GetProperty().SetLineWidth(2.0) + renderer.AddActor(actor) + + +def get_track(rec, pidx): + """Extract valid steps from record array.""" + steps = rec[pidx] + flagmasks = steps.view(np.uint32)[:, 3, 3] + n = int(np.sum(flagmasks != 0)) + return steps[:n, 0, :3], n + + +def add_track(renderer, track, n, color, radius=1.5, opacity=0.8): + """Add a polyline tube for a photon track.""" + if n < 2: + return + points = vtk.vtkPoints() + for i in range(n): + points.InsertNextPoint(float(track[i, 0]), float(track[i, 1]), float(track[i, 2])) + line = vtk.vtkPolyLine() + line.GetPointIds().SetNumberOfIds(n) + for i in range(n): + line.GetPointIds().SetId(i, i) + lines = vtk.vtkCellArray() + lines.InsertNextCell(line) + pd = vtk.vtkPolyData() + pd.SetPoints(points) + pd.SetLines(lines) + tube = vtk.vtkTubeFilter() + tube.SetInputData(pd) + tube.SetRadius(radius) + tube.SetNumberOfSides(8) + tube.Update() + mapper = vtk.vtkPolyDataMapper() + mapper.SetInputData(tube.GetOutput()) + actor = vtk.vtkActor() + actor.SetMapper(mapper) + actor.GetProperty().SetColor(*color) + actor.GetProperty().SetOpacity(opacity) + renderer.AddActor(actor) + + +def add_sphere(renderer, pos, color, radius=4): + s = vtk.vtkSphereSource() + s.SetCenter(float(pos[0]), float(pos[1]), float(pos[2])) + s.SetRadius(radius) + s.SetThetaResolution(16) + s.SetPhiResolution(16) + m = vtk.vtkPolyDataMapper() + m.SetInputConnection(s.GetOutputPort()) + a = vtk.vtkActor() + a.SetMapper(m) + a.GetProperty().SetColor(*color) + renderer.AddActor(a) + + +def add_label(renderer, text, x, y, color, size): + txt = vtk.vtkTextActor() + txt.SetInput(text) + txt.GetPositionCoordinate().SetCoordinateSystemToNormalizedViewport() + txt.GetPositionCoordinate().SetValue(x, y) + txt.GetTextProperty().SetColor(*color) + txt.GetTextProperty().SetFontSize(size) + txt.GetTextProperty().SetBold(True) + txt.GetTextProperty().SetFontFamilyToCourier() + renderer.AddViewProp(txt) + + +def main(): + if len(sys.argv) < 6: + print(f"Usage: {sys.argv[0]} [num_tracks] [output.png]") + sys.exit(1) + + gpu_rec = np.load(sys.argv[1]) + g4_rec = np.load(sys.argv[2]) + gpu_phot = np.load(sys.argv[3]) + g4_phot = np.load(sys.argv[4]) + gdml_path = sys.argv[5] + num_tracks = int(sys.argv[6]) if len(sys.argv) > 6 else 20 + output = sys.argv[7] if len(sys.argv) > 7 else "divergent_tracks.png" + + # Find divergent photons + gpu_flags = gpu_phot.view(np.uint32).reshape(-1, 4, 4)[:, 3, 0] & 0xFFFF + g4_flags = g4_phot.view(np.uint32).reshape(-1, 4, 4)[:, 3, 0] & 0xFFFF + divergent = np.where(gpu_flags != g4_flags)[0][:num_tracks] + print(f"Plotting {len(divergent)} divergent photons: {divergent}") + + # Parse geometry + solids, physvols, volumes = parse_gdml_boxes(gdml_path) + + # VTK setup + renderer = vtk.vtkRenderer() + renderer.SetBackground(0.95, 0.95, 0.95) + renWin = vtk.vtkRenderWindow() + renWin.SetOffScreenRendering(True) + renWin.SetSize(1800, 1200) + renWin.AddRenderer(renderer) + + add_geometry(renderer, solids, physvols, volumes) + + # Start marker + add_sphere(renderer, (0, 0, 0), (0, 0.8, 0), 6) + + # Draw tracks + for pidx in divergent: + gpu_track, gpu_n = get_track(gpu_rec, pidx) + g4_track, g4_n = get_track(g4_rec, pidx) + + add_track(renderer, gpu_track, gpu_n, (1, 0, 0), 1.2, 0.7) + add_track(renderer, g4_track, g4_n, (0.1, 0.3, 1), 1.2, 0.7) + + if gpu_n >= 2: + add_sphere(renderer, gpu_track[-1], (1, 0, 0), 4) + if g4_n >= 2: + add_sphere(renderer, g4_track[-1], (0.1, 0.3, 1), 4) + + # Labels + add_label(renderer, f"{len(divergent)} divergent photons: GPU vs G4", 0.02, 0.94, (0, 0, 0), 18) + add_label(renderer, "Red = GPU tracks", 0.02, 0.89, (0.8, 0, 0), 15) + add_label(renderer, "Blue = G4 tracks", 0.02, 0.84, (0, 0.2, 0.8), 15) + add_label(renderer, "Green = Start (0,0,0)", 0.02, 0.79, (0, 0.5, 0), 15) + + # Camera + camera = renderer.GetActiveCamera() + camera.SetPosition(600, -600, 500) + camera.SetFocalPoint(0, 0, 150) + camera.SetViewUp(0, 0, 1) + camera.SetViewAngle(45) + renderer.ResetCameraClippingRange() + + # Render and save + renWin.Render() + w2i = vtk.vtkWindowToImageFilter() + w2i.SetInput(renWin) + w2i.Update() + writer = vtk.vtkPNGWriter() + writer.SetFileName(output) + writer.SetInputConnection(w2i.GetOutputPort()) + writer.Write() + print(f"Saved {output}") + + +if __name__ == "__main__": + main() From 6b6f9f51955bd80e5b9bff15085d2f5b592e702f Mon Sep 17 00:00:00 2001 From: ggalgoczi Date: Fri, 10 Apr 2026 01:11:57 +0000 Subject: [PATCH 4/4] Add apex.gdml and sphere_leak_simple.gdml test geometries MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit apex.gdml: scintillator detector with WLS plates, SiPMs, and Vikuiti reflective foil wrapping (500x500x6mm flat detector). sphere_leak_simple.gdml: nested sphere mirror leak test with dielectric_metal skinsurface, no RINDEX on mirror material — verified zero-leak configuration for both GPU and G4. --- examples/AlignGpuCpuTracks/apex.gdml | 379 ++++++++++++++++++ .../AlignGpuCpuTracks/sphere_leak_simple.gdml | 96 +++++ 2 files changed, 475 insertions(+) create mode 100644 examples/AlignGpuCpuTracks/apex.gdml create mode 100644 examples/AlignGpuCpuTracks/sphere_leak_simple.gdml diff --git a/examples/AlignGpuCpuTracks/apex.gdml b/examples/AlignGpuCpuTracks/apex.gdml new file mode 100644 index 000000000..49d36656b --- /dev/null +++ b/examples/AlignGpuCpuTracks/apex.gdml @@ -0,0 +1,379 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/AlignGpuCpuTracks/sphere_leak_simple.gdml b/examples/AlignGpuCpuTracks/sphere_leak_simple.gdml new file mode 100644 index 000000000..f75a535f6 --- /dev/null +++ b/examples/AlignGpuCpuTracks/sphere_leak_simple.gdml @@ -0,0 +1,96 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +