From 8ce578f0e667661ab86c7228a67cd585c6740c69 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 1 May 2026 17:47:59 +0000 Subject: [PATCH 1/2] feat(examples): add GPURaytraceQuasi using G4 11.4 native offload classes New standalone example demonstrating optical-photon offload via the G4 11.4 G4QuasiCerenkov / G4QuasiScintillation processes. Mirrors the structure of examples/async_gpu_std: standalone CMakeLists with find_package(eic-opticks), bare includes, manual argv parsing, glm find_package workaround for the eic-opticksConfig.cmake gap. Quasi process registration uses an inline QuasiOpticalPhysics constructor because G4OpticalPhysics::ConstructProcess in v11.4.1 does not auto-instantiate the Quasi variants based on G4OpticalParameters offload flags. The constructor activates after G4OpticalPhysics has run, with the legacy Cerenkov / Scintillation processes disabled via SetProcessActivation(false). The SteppingAction matches both "QuasiScintillation" and the upstream typo "QausiScintillation" for forward compatibility. Bit-identical output to the legacy GPURaytrace example when both run with stackPhotons=false (verified at 100-event and 50M-photon scales). No modifications to eic-opticks core source. --- examples/GPURaytraceQuasi/CMakeLists.txt | 24 + .../GPURaytraceQuasi/GPURaytraceQuasi.cpp | 234 +++++++ examples/GPURaytraceQuasi/GPURaytraceQuasi.h | 626 ++++++++++++++++++ .../GPURaytraceQuasi/macros/run_parity.mac | 5 + .../GPURaytraceQuasi/macros/run_perf_50M.mac | 12 + .../GPURaytraceQuasi/macros/run_perf_fair.mac | 13 + 6 files changed, 914 insertions(+) create mode 100644 examples/GPURaytraceQuasi/CMakeLists.txt create mode 100644 examples/GPURaytraceQuasi/GPURaytraceQuasi.cpp create mode 100644 examples/GPURaytraceQuasi/GPURaytraceQuasi.h create mode 100644 examples/GPURaytraceQuasi/macros/run_parity.mac create mode 100644 examples/GPURaytraceQuasi/macros/run_perf_50M.mac create mode 100644 examples/GPURaytraceQuasi/macros/run_perf_fair.mac diff --git a/examples/GPURaytraceQuasi/CMakeLists.txt b/examples/GPURaytraceQuasi/CMakeLists.txt new file mode 100644 index 000000000..8441b1ffa --- /dev/null +++ b/examples/GPURaytraceQuasi/CMakeLists.txt @@ -0,0 +1,24 @@ +cmake_minimum_required(VERSION 3.20) + +project(GPURaytraceQuasi LANGUAGES CXX) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) + +# Order matters: pull in glm before eic-opticks because the installed +# eic-opticks::SysRap target's link interface references glm::glm but the +# package's Config.cmake does not declare it as a dependency. +find_package(glm REQUIRED) +find_package(eic-opticks REQUIRED) +find_package(Geant4 REQUIRED ui_all vis_all) + +add_executable(GPURaytraceQuasi GPURaytraceQuasi.cpp GPURaytraceQuasi.h) +target_link_libraries(GPURaytraceQuasi + eic-opticks::G4CX + eic-opticks::SysRap + eic-opticks::U4 + ${Geant4_LIBRARIES} +) + +install(TARGETS GPURaytraceQuasi) diff --git a/examples/GPURaytraceQuasi/GPURaytraceQuasi.cpp b/examples/GPURaytraceQuasi/GPURaytraceQuasi.cpp new file mode 100644 index 000000000..72123e683 --- /dev/null +++ b/examples/GPURaytraceQuasi/GPURaytraceQuasi.cpp @@ -0,0 +1,234 @@ +#include +#include +#include +#include + +#include "FTFP_BERT.hh" +#include "G4LossTableManager.hh" +#include "G4OpticalParameters.hh" +#include "G4OpticalPhysics.hh" +#include "G4ParticleTable.hh" +#include "G4ProcessManager.hh" +#include "G4QuasiCerenkov.hh" +#include "G4QuasiScintillation.hh" +#include "G4VModularPhysicsList.hh" +#include "G4VPhysicsConstructor.hh" + +#include "G4UIExecutive.hh" +#include "G4UImanager.hh" +#include "G4VisExecutive.hh" + +#include "sysrap/OPTICKS_LOG.hh" + +#include "GPURaytraceQuasi.h" + +// Companion physics constructor: replaces G4Cerenkov / G4Scintillation with their +// G4 11.4 Quasi variants. Registered AFTER G4OpticalPhysics; G4OpticalPhysics is +// configured to skip legacy Cerenkov / Scintillation registration via +// SetProcessActivation("Cerenkov", false) / SetProcessActivation("Scintillation", false). +// +// Required because G4OpticalPhysics::ConstructProcess in v11.4.1 hardcodes legacy +// G4Cerenkov / G4Scintillation and does not honour the SetCerenkovOffloadPhotons / +// SetScintOffloadPhotons flags. Verified against G4 v11.4.1 release source at +// source/physics_lists/constructors/electromagnetic/src/G4OpticalPhysics.cc — no +// branch on the offload parameter exists. +class QuasiOpticalPhysics : public G4VPhysicsConstructor +{ + public: + QuasiOpticalPhysics(const G4String &name = "QuasiOptical") : G4VPhysicsConstructor(name) {} + ~QuasiOpticalPhysics() override = default; + + void ConstructParticle() override + { + // G4OpticalPhysics already calls G4QuasiOpticalPhoton::QuasiOpticalPhotonDefinition() + // in its ConstructParticle(). Nothing to do. + } + + void ConstructProcess() override + { + auto *params = G4OpticalParameters::Instance(); + + G4QuasiCerenkov *theCerenkov = nullptr; + if (params->GetCerenkovOffloadPhotons()) + { + theCerenkov = new G4QuasiCerenkov(); + theCerenkov->SetOffloadPhotons(true); + theCerenkov->SetStackPhotons(params->GetCerenkovStackPhotons()); + } + + G4QuasiScintillation *theScint = nullptr; + if (params->GetScintOffloadPhotons()) + { + theScint = new G4QuasiScintillation(); + theScint->SetOffloadPhotons(true); + theScint->SetStackPhotons(params->GetScintStackPhotons()); + G4EmSaturation *emSat = G4LossTableManager::Instance()->EmSaturation(); + theScint->AddSaturation(emSat); + } + + auto *iter = GetParticleIterator(); + iter->reset(); + while ((*iter)()) + { + auto *particle = iter->value(); + if (particle->IsShortLived()) + continue; + auto *pManager = particle->GetProcessManager(); + if (!pManager) + continue; + + if (theCerenkov && theCerenkov->IsApplicable(*particle)) + { + pManager->AddDiscreteProcess(theCerenkov); + } + if (theScint && theScint->IsApplicable(*particle)) + { + pManager->AddProcess(theScint); + pManager->SetProcessOrderingToLast(theScint, idxAtRest); + pManager->SetProcessOrderingToLast(theScint, idxPostStep); + } + } + } +}; + +#include "G4RunManager.hh" +#include "G4RunManagerFactory.hh" +#include "G4VUserActionInitialization.hh" + +using namespace std; + +struct ActionInitialization : public G4VUserActionInitialization +{ + private: + G4App *fG4App; // Store the pointer to G4App + + public: + // Note the signature: now we take a pointer to the G4App itself + ActionInitialization(G4App *app) : G4VUserActionInitialization(), fG4App(app) + { + } + + virtual void BuildForMaster() const override + { + SetUserAction(fG4App->run_act_); + } + + virtual void Build() const override + { + SetUserAction(fG4App->prim_gen_); + SetUserAction(fG4App->run_act_); + SetUserAction(fG4App->event_act_); + SetUserAction(fG4App->tracking_); + SetUserAction(fG4App->stepping_); + } +}; + +static void usage(const char *prog) +{ + std::cerr << "Usage: " << prog + << " [options]\n" + " -g, --gdml PATH GDML file (default: geom.gdml)\n" + " -m, --macro PATH Geant4 macro (default: run.mac)\n" + " -s, --seed N random seed (default: time())\n" + " -i, --interactive open interactive viewer\n" + " -h, --help show this message\n"; +} + +int main(int argc, char **argv) +{ + OPTICKS_LOG(argc, argv); + + string gdml_file = "geom.gdml"; + string macro_name = "run.mac"; + long seed = static_cast(std::time(nullptr)); + bool interactive = false; + + for (int i = 1; i < argc; i++) + { + string a = argv[i]; + auto next = [&](const char *flag) -> const char * { + if (i + 1 >= argc) + { + std::cerr << flag << " requires an argument\n"; + std::exit(EXIT_FAILURE); + } + return argv[++i]; + }; + + if (a == "-g" || a == "--gdml") + gdml_file = next("--gdml"); + else if (a == "-m" || a == "--macro") + macro_name = next("--macro"); + else if (a == "-s" || a == "--seed") + seed = std::atol(next("--seed")); + else if (a == "-i" || a == "--interactive") + interactive = true; + else if (a == "-h" || a == "--help") + { + usage(argv[0]); + return EXIT_SUCCESS; + } + else + { + std::cerr << "unknown option: " << a << "\n"; + usage(argv[0]); + return EXIT_FAILURE; + } + } + + CLHEP::HepRandom::setTheSeed(seed); + G4cout << "Random seed set to: " << seed << G4endl; + + // Configure Geant4 optical photon offload (G4 11.4 G4QuasiCerenkov / G4QuasiScintillation). + // + // G4OpticalPhysics::ConstructProcess in v11.4.1 (and current master) hardcodes legacy + // G4Cerenkov / G4Scintillation and does not branch on the offload flags. So: + // 1) Tell G4OpticalPhysics to skip Cerenkov / Scintillation entirely. + // 2) Register a small companion physics constructor that installs the Quasi variants. + // + // StackPhotons=false: no real G4OpticalPhoton secondaries; Opticks does the propagation. + // OffloadPhotons=true: read by our QuasiOpticalPhysics constructor to gate registration. + auto *optParams = G4OpticalParameters::Instance(); + optParams->SetProcessActivation("Cerenkov", false); + optParams->SetProcessActivation("Scintillation", false); + optParams->SetCerenkovOffloadPhotons(true); + optParams->SetCerenkovStackPhotons(false); + optParams->SetScintOffloadPhotons(true); + optParams->SetScintStackPhotons(false); + + // The physics list must be instantiated before other user actions + G4VModularPhysicsList *physics = new FTFP_BERT; + physics->RegisterPhysics(new G4OpticalPhysics); + physics->RegisterPhysics(new QuasiOpticalPhysics); + + auto *run_mgr = G4RunManagerFactory::CreateRunManager(); + run_mgr->SetUserInitialization(physics); + + G4App *g4app = new G4App(gdml_file); + + ActionInitialization *actionInit = new ActionInitialization(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; +} diff --git a/examples/GPURaytraceQuasi/GPURaytraceQuasi.h b/examples/GPURaytraceQuasi/GPURaytraceQuasi.h new file mode 100644 index 000000000..c11ae0e41 --- /dev/null +++ b/examples/GPURaytraceQuasi/GPURaytraceQuasi.h @@ -0,0 +1,626 @@ +#include +#include +#include + +#include "G4AutoLock.hh" +#include "G4BooleanSolid.hh" +#include "G4Cerenkov.hh" +#include "G4Electron.hh" +#include "G4Event.hh" +#include "G4GDMLParser.hh" +#include "G4LogicalVolumeStore.hh" +#include "G4OpBoundaryProcess.hh" +#include "G4OpticalPhoton.hh" +#include "G4PhysicalConstants.hh" +#include "G4PrimaryParticle.hh" +#include "G4PrimaryVertex.hh" +#include "G4QuasiCerenkov.hh" +#include "G4QuasiOpticalPhoton.hh" +#include "G4QuasiScintillation.hh" +#include "G4RunManager.hh" +#include "G4RunManagerFactory.hh" +#include "G4SDManager.hh" +#include "G4Scintillation.hh" +#include "G4SubtractionSolid.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 "G4VProcess.hh" +#include "G4VUserDetectorConstruction.hh" +#include "G4VUserPrimaryGeneratorAction.hh" + +#include "g4cx/G4CXOpticks.hh" +#include "sysrap/NP.hh" +#include "sysrap/SEvt.hh" +#include "sysrap/STrackInfo.h" +#include "sysrap/spho.h" +#include "sysrap/sphoton.h" +#include "u4/U4.hh" +#include "u4/U4Random.hh" +#include "u4/U4StepPoint.hh" +#include "u4/U4Touchable.h" +#include "u4/U4Track.h" + +namespace +{ +G4Mutex genstep_mutex = G4MUTEX_INITIALIZER; +} + +bool IsSubtractionSolid(G4VSolid *solid) +{ + if (!solid) + return false; + + // Check if the solid is directly a G4SubtractionSolid + if (dynamic_cast(solid)) + return true; + + // If the solid is a Boolean solid, check its constituent solids + G4BooleanSolid *booleanSolid = dynamic_cast(solid); + if (booleanSolid) + { + G4VSolid *solidA = booleanSolid->GetConstituentSolid(0); + G4VSolid *solidB = booleanSolid->GetConstituentSolid(1); + + // Recursively check the constituent solids + if (IsSubtractionSolid(solidA) || IsSubtractionSolid(solidB)) + return true; + } + + // For other solid types, return false + return false; +} + +std::string str_tolower(std::string s) +{ + std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c) { return std::tolower(c); }); + return s; +} + +struct PhotonHit : public G4VHit +{ + PhotonHit() = default; + + PhotonHit(unsigned id, G4double energy, G4double time, G4ThreeVector position, G4ThreeVector direction, + G4ThreeVector polarization) + : fid(id), fenergy(energy), ftime(time), fposition(position), fdirection(direction), fpolarization(polarization) + { + } + + // Copy constructor + PhotonHit(const PhotonHit &right) + : G4VHit(right), fid(right.fid), fenergy(right.fenergy), ftime(right.ftime), fposition(right.fposition), + fdirection(right.fdirection), fpolarization(right.fpolarization) + { + } + + // Assignment operator + const PhotonHit &operator=(const PhotonHit &right) + { + if (this != &right) + { + G4VHit::operator=(right); + fid = right.fid; + fenergy = right.fenergy; + ftime = right.ftime; + fposition = right.fposition; + fdirection = right.fdirection; + fpolarization = right.fpolarization; + } + return *this; + } + + // Equality operator + G4bool operator==(const PhotonHit &right) const + { + return (this == &right); + } + + // Print method + void Print() override + { + G4cout << "Detector id: " << fid << " energy: " << fenergy << " nm" << " time: " << ftime << " ns" + << " position: " << fposition << " direction: " << fdirection << " polarization: " << fpolarization + << G4endl; + } + + // Member variables + G4int fid{0}; + G4double fenergy{0}; + G4double ftime{0}; + G4ThreeVector fposition{0, 0, 0}; + G4ThreeVector fdirection{0, 0, 0}; + G4ThreeVector fpolarization{0, 0, 0}; +}; + +using PhotonHitsCollection = G4THitsCollection; + +struct PhotonSD : public G4VSensitiveDetector +{ + PhotonSD(G4String name) : G4VSensitiveDetector(name), fHCID(-1) + { + G4String HCname = name + "_HC"; + collectionName.insert(HCname); + G4cout << collectionName.size() << " PhotonSD name: " << name << " collection Name: " << HCname << G4endl; + } + + void Initialize(G4HCofThisEvent *hce) override + { + fPhotonHitsCollection = new PhotonHitsCollection(SensitiveDetectorName, collectionName[0]); + if (fHCID < 0) + { + // G4cout << "PhotonSD::Initialize: " << SensitiveDetectorName << " " << collectionName[0] << G4endl; + fHCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); + } + hce->AddHitsCollection(fHCID, fPhotonHitsCollection); + } + + G4bool ProcessHits(G4Step *aStep, G4TouchableHistory *) override + { + G4Track *theTrack = aStep->GetTrack(); + if (theTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) + return false; + + G4double theEnergy = theTrack->GetTotalEnergy() / CLHEP::eV; + + // Create a new hit (CopyNr is set to 0 as DetectorID is omitted) + PhotonHit *newHit = new PhotonHit( + 0, // CopyNr set to 0 + theEnergy, theTrack->GetGlobalTime(), aStep->GetPostStepPoint()->GetPosition(), + aStep->GetPostStepPoint()->GetMomentumDirection(), aStep->GetPostStepPoint()->GetPolarization()); + + fPhotonHitsCollection->insert(newHit); + theTrack->SetTrackStatus(fStopAndKill); + return true; + } + + void EndOfEvent(G4HCofThisEvent *) override + { + G4int NbHits = fPhotonHitsCollection->entries(); + G4cout << "PhotonSD::EndOfEvent Number of PhotonHits: " << NbHits << G4endl; + } + + void AddOpticksHits() + { + SEvt *sev = SEvt::Get_EGPU(); + unsigned int num_hits = sev->GetNumHit(0); + + for (int idx = 0; idx < int(num_hits); idx++) + { + sphoton hit; + sev->getHit(hit, idx); + G4ThreeVector position = G4ThreeVector(hit.pos.x, hit.pos.y, hit.pos.z); + G4ThreeVector direction = G4ThreeVector(hit.mom.x, hit.mom.y, hit.mom.z); + G4ThreeVector polarization = G4ThreeVector(hit.pol.x, hit.pol.y, hit.pol.z); + int theCreationProcessid; + if (OpticksPhoton::HasCerenkovFlag(hit.flagmask)) + { + theCreationProcessid = 0; + } + else if (OpticksPhoton::HasScintillationFlag(hit.flagmask)) + { + theCreationProcessid = 1; + } + else + { + theCreationProcessid = -1; + } + std::cout << hit.wavelength << " " << position << " " << direction << " " << polarization << std::endl; + + PhotonHit *newHit = new PhotonHit(0, hit.wavelength, hit.time, position, direction, polarization); + fPhotonHitsCollection->insert(newHit); + } + } + + private: + PhotonHitsCollection *fPhotonHitsCollection{nullptr}; + G4int fHCID; +}; + +struct DetectorConstruction : G4VUserDetectorConstruction +{ + DetectorConstruction(std::filesystem::path gdml_file) : gdml_file_(gdml_file) + { + } + + G4VPhysicalVolume *Construct() override + { + parser_.Read(gdml_file_.string(), false); + G4VPhysicalVolume *world = parser_.GetWorldVolume(); + + G4CXOpticks::SetGeometry(world); + G4LogicalVolumeStore *lvStore = G4LogicalVolumeStore::GetInstance(); + + static G4VisAttributes invisibleVisAttr(false); + + // Check if the store is not empty + if (lvStore && !lvStore->empty()) + { + // Iterate over all logical volumes in the store + for (auto &logicalVolume : *lvStore) + { + G4VSolid *solid = logicalVolume->GetSolid(); + + // Check if the solid uses subtraction + if (IsSubtractionSolid(solid)) + { + // Assign the invisible visual attributes to the logical volume + logicalVolume->SetVisAttributes(&invisibleVisAttr); + + // Optionally, print out the name of the logical volume + G4cout << "Hiding logical volume: " << logicalVolume->GetName() << G4endl; + } + } + } + + return world; + } + + void ConstructSDandField() override + { + G4cout << "ConstructSDandField is called." << G4endl; + G4SDManager *SDman = G4SDManager::GetSDMpointer(); + + const G4GDMLAuxMapType *auxmap = parser_.GetAuxMap(); + for (auto const &[logVol, listType] : *auxmap) + { + for (auto const &auxtype : listType) + { + if (auxtype.type == "SensDet") + { + G4cout << "Attaching sensitive detector to logical volume: " << logVol->GetName() << G4endl; + G4String name = logVol->GetName() + "_PhotonDetector"; + PhotonSD *aPhotonSD = new PhotonSD(name); + SDman->AddNewDetector(aPhotonSD); + logVol->SetSensitiveDetector(aPhotonSD); + } + } + } + } + + private: + std::filesystem::path gdml_file_; + G4GDMLParser parser_; +}; + +struct PrimaryGenerator : G4VUserPrimaryGeneratorAction +{ + SEvt *sev; + + PrimaryGenerator(SEvt *sev) : sev(sev) + { + } + + void GeneratePrimaries(G4Event *event) override + { + G4ThreeVector position_mm(0.0 * m, 0.0 * m, 0.0 * m); + G4double time_ns = 0; + G4ThreeVector direction(0, 0.2, 0.8); + G4double wavelength_nm = 0.1; + + G4PrimaryVertex *vertex = new G4PrimaryVertex(position_mm, time_ns); + G4PrimaryParticle *particle = new G4PrimaryParticle(G4Electron::Definition()); + particle->SetKineticEnergy(5 * GeV); + particle->SetMomentumDirection(direction); + vertex->SetPrimary(particle); + event->AddPrimaryVertex(vertex); + } +}; + +struct EventAction : G4UserEventAction +{ + SEvt *sev; + G4int fTotalG4Hits{0}; + + EventAction(SEvt *sev) : sev(sev) + { + } + + void BeginOfEventAction(const G4Event *event) override + { + } + + void EndOfEventAction(const G4Event *event) override + { + G4HCofThisEvent *hce = event->GetHCofThisEvent(); + if (hce) + { + for (G4int i = 0; i < hce->GetNumberOfCollections(); i++) + { + G4VHitsCollection *hc = hce->GetHC(i); + if (hc) + { + fTotalG4Hits += hc->GetSize(); + } + } + } + } + + G4int GetTotalG4Hits() const + { + return fTotalG4Hits; + } +}; + +struct RunAction : G4UserRunAction +{ + EventAction *fEventAction; + + RunAction(EventAction *eventAction) : fEventAction(eventAction) + { + } + + void BeginOfRunAction(const G4Run *run) override + { + } + + void EndOfRunAction(const G4Run *run) override + { + if (G4Threading::IsMasterThread()) + { + 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(); + // Compute duration + std::chrono::duration elapsed = end - start; + std::cout << "Simulation time: " << elapsed.count() << " seconds" << std::endl; + + // unsigned int num_hits = SEvt::GetNumHit(EGPU); + SEvt *sev = SEvt::Get_EGPU(); + unsigned int num_hits = sev->GetNumHit(0); + + std::cout << "Opticks: NumCollected: " << sev->GetNumGenstepFromGenstep(0) << std::endl; + std::cout << "Opticks: NumCollected: " << sev->GetNumPhotonCollected(0) << std::endl; + std::cout << "Opticks: NumHits: " << num_hits << std::endl; + std::cout << "Geant4: NumHits: " << fEventAction->GetTotalG4Hits() << std::endl; + std::ofstream outFile("opticks_hits_output.txt"); + if (!outFile.is_open()) + { + std::cerr << "Error opening output file!" << std::endl; + return; + } + + for (int idx = 0; idx < int(num_hits); idx++) + { + sphoton hit; + sev->getHit(hit, idx); + G4ThreeVector position = G4ThreeVector(hit.pos.x, hit.pos.y, hit.pos.z); + G4ThreeVector direction = G4ThreeVector(hit.mom.x, hit.mom.y, hit.mom.z); + G4ThreeVector polarization = G4ThreeVector(hit.pol.x, hit.pol.y, hit.pol.z); + int theCreationProcessid; + if (OpticksPhoton::HasCerenkovFlag(hit.flagmask)) + { + theCreationProcessid = 0; + } + else if (OpticksPhoton::HasScintillationFlag(hit.flagmask)) + { + theCreationProcessid = 1; + } + else + { + theCreationProcessid = -1; + } + // std::cout << "Adding hit from Opticks:" << hit.wavelength << " " << position << " " << direction + // << " + // " + // << polarization << std::endl; + outFile << hit.time << " " << hit.wavelength << " " << "(" << position.x() << ", " << position.y() + << ", " << position.z() << ") " << "(" << direction.x() << ", " << direction.y() << ", " + << direction.z() << ") " << "(" << polarization.x() << ", " << polarization.y() << ", " + << polarization.z() << ") " << "CreationProcessID=" << theCreationProcessid << std::endl; + } + + outFile.close(); + } + } +}; + +struct SteppingAction : G4UserSteppingAction +{ + SEvt *sev; + + SteppingAction(SEvt *sev) : sev(sev) + { + } + + void UserSteppingAction(const G4Step *aStep) + { + G4Track *aTrack; + G4int fNumPhotons = 0; + + G4StepPoint *preStep = aStep->GetPostStepPoint(); + G4VPhysicalVolume *volume = preStep->GetPhysicalVolume(); + + if (aStep->GetTrack()->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) + { + // Kill if step count exceeds 10000 to avoid reflection forever + if (aStep->GetTrack()->GetCurrentStepNumber() > 10000) + { + aStep->GetTrack()->SetTrackStatus(fStopAndKill); + } + } + + // G4QuasiCerenkov / G4QuasiScintillation create a G4QuasiOpticalPhoton secondary + // per parent step (carrying the burst metadata as aux info). We capture the + // genstep at the parent's step in the QuasiCerenkov / QuasiScintillation branches + // below, so the secondary itself has no further purpose — kill it to avoid + // tracking through a particle with no registered processes. + if (aStep->GetTrack()->GetDefinition() == G4QuasiOpticalPhoton::Definition()) + { + aStep->GetTrack()->SetTrackStatus(fStopAndKill); + return; + } + + if (volume && volume->GetName() == "MirrorPyramid") + { + aTrack = aStep->GetTrack(); + if (aTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) + { + aTrack->SetTrackStatus(fStopAndKill); + } + } + + G4SteppingManager *fpSteppingManager = + G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager(); + G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus(); + + if (stepStatus != fAtRestDoItProc) + { + G4ProcessVector *procPost = fpSteppingManager->GetfPostStepDoItVector(); + size_t MAXofPostStepLoops = fpSteppingManager->GetMAXofPostStepLoops(); + for (size_t i3 = 0; i3 < MAXofPostStepLoops; i3++) + { + const G4String &procName = (*procPost)[i3]->GetProcessName(); + // Diagnostic (Step 4 fold-in): once per worker, dump the optical processes + // we see to verify G4OpticalPhysics installed the Quasi variants and not + // the legacy ones. Removed in Step 8. + static thread_local bool quasi_dump_done = false; + if (!quasi_dump_done && + (procName == "Cerenkov" || procName == "QuasiCerenkov" || + procName == "Scintillation" || procName == "QuasiScintillation" || + procName == "QausiScintillation")) + { + G4cout << "[GPURaytraceQuasi] post-step process found: " << procName << G4endl; + if (procName == "Cerenkov" || procName == "Scintillation") + { + G4cout << "[GPURaytraceQuasi] WARNING: legacy " << procName + << " process is registered — offload swap may not have taken effect." + << G4endl; + } + } + + if (procName == "QuasiCerenkov") + { + aTrack = aStep->GetTrack(); + const G4DynamicParticle *aParticle = aTrack->GetDynamicParticle(); + G4double charge = aParticle->GetDefinition()->GetPDGCharge(); + const G4Material *aMaterial = aTrack->GetMaterial(); + G4MaterialPropertiesTable *MPT = aMaterial->GetMaterialPropertiesTable(); + + G4MaterialPropertyVector *Rindex = MPT->GetProperty(kRINDEX); + if (!Rindex || Rindex->GetVectorLength() == 0) + { + G4cout << "WARNING: Material has no valid RINDEX data. Skipping Cerenkov calculation." + << G4endl; + return; + } + + G4QuasiCerenkov *proc = (G4QuasiCerenkov *)(*procPost)[i3]; + fNumPhotons = proc->GetNumPhotons(); + + G4AutoLock lock(&genstep_mutex); // <-- Mutex is locked here + + if (fNumPhotons > 0) + { + G4double Pmin = Rindex->Energy(0); + G4double Pmax = Rindex->GetMaxEnergy(); + G4double nMax = Rindex->GetMaxValue(); + G4double beta1 = aStep->GetPreStepPoint()->GetBeta(); + G4double beta2 = aStep->GetPostStepPoint()->GetBeta(); + G4double beta = (beta1 + beta2) * 0.5; + G4double BetaInverse = 1. / beta; + G4double maxCos = BetaInverse / nMax; + G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos); + G4double MeanNumberOfPhotons1 = + proc->GetAverageNumberOfPhotons(charge, beta1, aMaterial, Rindex); + G4double MeanNumberOfPhotons2 = + proc->GetAverageNumberOfPhotons(charge, beta2, aMaterial, Rindex); + + U4::CollectGenstep_G4Cerenkov_modified(aTrack, aStep, fNumPhotons, BetaInverse, Pmin, Pmax, + maxCos, maxSin2, MeanNumberOfPhotons1, + MeanNumberOfPhotons2); + if (!quasi_dump_done) + { + G4cout << "[GPURaytraceQuasi] first QuasiCerenkov genstep captured: N=" + << fNumPhotons << G4endl; + quasi_dump_done = true; + } + } + } + if (procName == "QuasiScintillation" || procName == "QausiScintillation") + { + G4QuasiScintillation *proc1 = (G4QuasiScintillation *)(*procPost)[i3]; + fNumPhotons = proc1->GetNumPhotons(); + if (fNumPhotons > 0) + { + aTrack = aStep->GetTrack(); + const G4Material *aMaterial = aTrack->GetMaterial(); + G4MaterialPropertiesTable *MPT = aMaterial->GetMaterialPropertiesTable(); + if (!MPT || !MPT->ConstPropertyExists(kSCINTILLATIONTIMECONSTANT1)) + { + G4cout << "WARNING: Material has no valid SCINTILLATIONTIMECONSTANT1 data. Skipping " + "Scintillation calculation." + << G4endl; + return; + } + G4double SCINTILLATIONTIMECONSTANT1 = MPT->GetConstProperty(kSCINTILLATIONTIMECONSTANT1); + + U4::CollectGenstep_DsG4Scintillation_r4695(aTrack, aStep, fNumPhotons, 1, + SCINTILLATIONTIMECONSTANT1); + if (!quasi_dump_done) + { + G4cout << "[GPURaytraceQuasi] first QuasiScintillation genstep captured: N=" + << fNumPhotons << G4endl; + quasi_dump_done = true; + } + } + } + } + } + } +}; + +struct TrackingAction : G4UserTrackingAction +{ + const G4Track *transient_fSuspend_track = nullptr; + SEvt *sev; + + TrackingAction(SEvt *sev) : sev(sev) + { + } + + void PreUserTrackingAction_Optical_FabricateLabel(const G4Track *track) + { + } + + void PreUserTrackingAction(const G4Track *track) override + { + } + + void PostUserTrackingAction(const G4Track *track) override + { + } +}; + +struct G4App +{ + G4App(std::filesystem::path gdml_file) + : sev(SEvt::CreateOrReuse_EGPU()), det_cons_(new DetectorConstruction(gdml_file)), + prim_gen_(new PrimaryGenerator(sev)), event_act_(new EventAction(sev)), run_act_(new RunAction(event_act_)), + stepping_(new SteppingAction(sev)), + + tracking_(new TrackingAction(sev)) + { + } + + //~G4App(){ G4CXOpticks::Finalize();} + + // Create "global" event + SEvt *sev; + + G4VUserDetectorConstruction *det_cons_; + G4VUserPrimaryGeneratorAction *prim_gen_; + EventAction *event_act_; + RunAction *run_act_; + SteppingAction *stepping_; + TrackingAction *tracking_; +}; diff --git a/examples/GPURaytraceQuasi/macros/run_parity.mac b/examples/GPURaytraceQuasi/macros/run_parity.mac new file mode 100644 index 000000000..db2add916 --- /dev/null +++ b/examples/GPURaytraceQuasi/macros/run_parity.mac @@ -0,0 +1,5 @@ +# Parity-comparison macro: 100 electrons in one run to reduce per-event +# RNG-noise fluctuation in the OLD vs NEW (Quasi) hit-count comparison. +/run/verbose 1 +/run/initialize +/run/beamOn 100 diff --git a/examples/GPURaytraceQuasi/macros/run_perf_50M.mac b/examples/GPURaytraceQuasi/macros/run_perf_50M.mac new file mode 100644 index 000000000..7c289a5ab --- /dev/null +++ b/examples/GPURaytraceQuasi/macros/run_perf_50M.mac @@ -0,0 +1,12 @@ +# Performance test at ~50M photon scale. +# - single-thread (avoids the pre-existing eic-opticks MT crash on beamOn>1) +# - stackPhotons=false on both Cerenkov + Scintillation (Opticks-only mode) +# - beamOn 4000: ~14k Cerenkov photons per 1 GeV e- in raindrop +# → ~56M photons total for cerenkov-only test +# → ~80M+ photons for cerenkov+scintillation test +/run/numberOfThreads 1 +/process/optical/cerenkov/setStackPhotons false +/process/optical/scintillation/setStackPhotons false +/run/verbose 0 +/run/initialize +/run/beamOn 4000 diff --git a/examples/GPURaytraceQuasi/macros/run_perf_fair.mac b/examples/GPURaytraceQuasi/macros/run_perf_fair.mac new file mode 100644 index 000000000..da8ff62f7 --- /dev/null +++ b/examples/GPURaytraceQuasi/macros/run_perf_fair.mac @@ -0,0 +1,13 @@ +# Fair performance comparison macro: +# - single-thread (avoids the pre-existing eic-opticks MT crash on beamOn>1 +# with the scintillation geometry) +# - stackPhotons=false on both Cerenkov and Scintillation, so the legacy +# GPURaytrace and the new GPURaytraceQuasi both run pure offload mode +# (no real G4OpticalPhoton secondaries; Opticks-only propagation) +# - beamOn 100 so per-event statistical noise averages out +/run/numberOfThreads 1 +/process/optical/cerenkov/setStackPhotons false +/process/optical/scintillation/setStackPhotons false +/run/verbose 0 +/run/initialize +/run/beamOn 100 From c4f6f988918925489e9ce920783491bd9fcb1fa4 Mon Sep 17 00:00:00 2001 From: Gabor Galgoczi Date: Fri, 1 May 2026 17:54:49 +0000 Subject: [PATCH 2/2] style(GPURaytraceQuasi): apply repo clang-format style --- .../GPURaytraceQuasi/GPURaytraceQuasi.cpp | 9 ++- examples/GPURaytraceQuasi/GPURaytraceQuasi.h | 57 +++++++++++++------ 2 files changed, 46 insertions(+), 20 deletions(-) diff --git a/examples/GPURaytraceQuasi/GPURaytraceQuasi.cpp b/examples/GPURaytraceQuasi/GPURaytraceQuasi.cpp index 72123e683..ea091aa9c 100644 --- a/examples/GPURaytraceQuasi/GPURaytraceQuasi.cpp +++ b/examples/GPURaytraceQuasi/GPURaytraceQuasi.cpp @@ -35,7 +35,10 @@ class QuasiOpticalPhysics : public G4VPhysicsConstructor { public: - QuasiOpticalPhysics(const G4String &name = "QuasiOptical") : G4VPhysicsConstructor(name) {} + QuasiOpticalPhysics(const G4String &name = "QuasiOptical") : + G4VPhysicsConstructor(name) + { + } ~QuasiOpticalPhysics() override = default; void ConstructParticle() override @@ -104,7 +107,9 @@ struct ActionInitialization : public G4VUserActionInitialization public: // Note the signature: now we take a pointer to the G4App itself - ActionInitialization(G4App *app) : G4VUserActionInitialization(), fG4App(app) + ActionInitialization(G4App *app) : + G4VUserActionInitialization(), + fG4App(app) { } diff --git a/examples/GPURaytraceQuasi/GPURaytraceQuasi.h b/examples/GPURaytraceQuasi/GPURaytraceQuasi.h index c11ae0e41..26e8f2d15 100644 --- a/examples/GPURaytraceQuasi/GPURaytraceQuasi.h +++ b/examples/GPURaytraceQuasi/GPURaytraceQuasi.h @@ -88,15 +88,25 @@ struct PhotonHit : public G4VHit PhotonHit() = default; PhotonHit(unsigned id, G4double energy, G4double time, G4ThreeVector position, G4ThreeVector direction, - G4ThreeVector polarization) - : fid(id), fenergy(energy), ftime(time), fposition(position), fdirection(direction), fpolarization(polarization) + G4ThreeVector polarization) : + fid(id), + fenergy(energy), + ftime(time), + fposition(position), + fdirection(direction), + fpolarization(polarization) { } // Copy constructor - PhotonHit(const PhotonHit &right) - : G4VHit(right), fid(right.fid), fenergy(right.fenergy), ftime(right.ftime), fposition(right.fposition), - fdirection(right.fdirection), fpolarization(right.fpolarization) + PhotonHit(const PhotonHit &right) : + G4VHit(right), + fid(right.fid), + fenergy(right.fenergy), + ftime(right.ftime), + fposition(right.fposition), + fdirection(right.fdirection), + fpolarization(right.fpolarization) { } @@ -143,7 +153,9 @@ using PhotonHitsCollection = G4THitsCollection; struct PhotonSD : public G4VSensitiveDetector { - PhotonSD(G4String name) : G4VSensitiveDetector(name), fHCID(-1) + PhotonSD(G4String name) : + G4VSensitiveDetector(name), + fHCID(-1) { G4String HCname = name + "_HC"; collectionName.insert(HCname); @@ -225,7 +237,8 @@ struct PhotonSD : public G4VSensitiveDetector struct DetectorConstruction : G4VUserDetectorConstruction { - DetectorConstruction(std::filesystem::path gdml_file) : gdml_file_(gdml_file) + DetectorConstruction(std::filesystem::path gdml_file) : + gdml_file_(gdml_file) { } @@ -293,7 +306,8 @@ struct PrimaryGenerator : G4VUserPrimaryGeneratorAction { SEvt *sev; - PrimaryGenerator(SEvt *sev) : sev(sev) + PrimaryGenerator(SEvt *sev) : + sev(sev) { } @@ -318,7 +332,8 @@ struct EventAction : G4UserEventAction SEvt *sev; G4int fTotalG4Hits{0}; - EventAction(SEvt *sev) : sev(sev) + EventAction(SEvt *sev) : + sev(sev) { } @@ -352,7 +367,8 @@ struct RunAction : G4UserRunAction { EventAction *fEventAction; - RunAction(EventAction *eventAction) : fEventAction(eventAction) + RunAction(EventAction *eventAction) : + fEventAction(eventAction) { } @@ -428,7 +444,8 @@ struct SteppingAction : G4UserSteppingAction { SEvt *sev; - SteppingAction(SEvt *sev) : sev(sev) + SteppingAction(SEvt *sev) : + sev(sev) { } @@ -584,7 +601,8 @@ struct TrackingAction : G4UserTrackingAction const G4Track *transient_fSuspend_track = nullptr; SEvt *sev; - TrackingAction(SEvt *sev) : sev(sev) + TrackingAction(SEvt *sev) : + sev(sev) { } @@ -603,12 +621,15 @@ struct TrackingAction : G4UserTrackingAction struct G4App { - G4App(std::filesystem::path gdml_file) - : sev(SEvt::CreateOrReuse_EGPU()), det_cons_(new DetectorConstruction(gdml_file)), - prim_gen_(new PrimaryGenerator(sev)), event_act_(new EventAction(sev)), run_act_(new RunAction(event_act_)), - stepping_(new SteppingAction(sev)), - - tracking_(new TrackingAction(sev)) + G4App(std::filesystem::path gdml_file) : + sev(SEvt::CreateOrReuse_EGPU()), + det_cons_(new DetectorConstruction(gdml_file)), + prim_gen_(new PrimaryGenerator(sev)), + event_act_(new EventAction(sev)), + run_act_(new RunAction(event_act_)), + stepping_(new SteppingAction(sev)), + + tracking_(new TrackingAction(sev)) { }