diff --git a/CSG/csg_intersect_leaf_newcone.h b/CSG/csg_intersect_leaf_newcone.h index 576c39bff3..6c335b2438 100644 --- a/CSG/csg_intersect_leaf_newcone.h +++ b/CSG/csg_intersect_leaf_newcone.h @@ -119,5 +119,3 @@ void intersect_leaf_newcone( bool& valid_isect, float4& isect, const quad& q0, c isect.w = t_cand ; } } - - diff --git a/README.md b/README.md index ab657e6818..ae63508f38 100644 --- a/README.md +++ b/README.md @@ -125,3 +125,66 @@ Process](https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/ ... ``` + + +## User/developer defined inputs + +### Defining primary particles + +There are certain user defined inputs that the user/developer has to define. In +the ```src/simg4oxmt``` example that imports ```src/g4appmt.h``` we provide +a working example with a simple geometry. The User/developer has to change the +following details: **Number of primary particles** to simulate in a macro file +and the **number of G4 threads**. For example: + +``` +/run/numberOfThreads {threads} +/run/verbose 1 +/process/optical/cerenkov/setStackPhotons {flag} +/run/initialize +/run/beamOn 500 +``` + +Here setStackPhotons defines **whether G4 will propagate optical photons or +not**. In production Opticks (GPU) takes care of the optical photon propagation. +Additionally the user has to define the **starting position**, **momentum** etc +of the primary particles define in the **GeneratePrimaries** function in +``src/g4appmt.h```. The hits of the optical photons are returned in the +**EndOfRunAction** function. If more photons are simulated than can fit in the +GPU RAM the execution of a GPU call should be moved to **EndOfEventAction** +together with retriving the hits. + +### Loading in geometry into EIC-Opticks + +EIC-Opticks can import geometries with GDML format automatically. There are +about 10 primitives supported now, eg. G4Box. G4Trd or G4Trap are not supported +yet, we are working on them. ```src/simg4oxmt``` takes GDML files through +arguments, eg. ```src/simg4oxmt -g mygdml.gdml```. + +The GDML must define all optical properties of surfaces of materials including: +- Efficiency (used by EIC-Opticks to specify detection efficiency and assign + sensitive surfaces) +- Refractive index +- Group velocity +- Reflectivity +- Etc. + + +## Performance studies + +In order to quantify the speed-up achieved by EIC-Opticks compared to G4 we +provide a python code that runs the same G4 simulation with and without tracking +optical photons in G4. The difference of the runs will yield the time required +to simulate photons. Meanwhile the same photons are simulated on GPU with +EIC-Opticks and the simulation time is saved. + +``` +mkdir -p /tmp/out/dev +mkdir -p /tmp/out/rel + +docker build -t eic-opticks:perf-dev --target=develop +docker run --rm -t -v /tmp/out:/tmp/out eic-opticks:perf-dev run-performance -o /tmp/out/dev + +docker build -t eic-opticks:perf-rel --target=release +docker run --rm -t -v /tmp/out:/tmp/out eic-opticks:perf-rel run-performance -o /tmp/out/rel +``` diff --git a/optiphy/tools/run_performance.py b/optiphy/tools/run_performance.py new file mode 100755 index 0000000000..c73fb1862c --- /dev/null +++ b/optiphy/tools/run_performance.py @@ -0,0 +1,114 @@ +import argparse +import subprocess +import re +import os +from pathlib import Path + +run_mac_template = """ +/run/numberOfThreads {threads} +/run/verbose 1 +/process/optical/cerenkov/setStackPhotons {flag} +/run/initialize +/run/beamOn 500 +""" + +os.environ["OPTICKS_EVENT_MODE"] = "Minimal" + +def get_opticks_home(): + """Get OPTICKS_HOME from environment, warn if not set.""" + opticks_home = os.environ.get("OPTICKS_HOME") + if opticks_home is None: + print("Warning: $OPTICKS_HOME is not defined, so this script should be called from the eic-opticks directory.") + return Path(".") + return Path(opticks_home) + +def parse_real_time(time_str): + # Parses 'real\t0m41.149s' to seconds + match = re.search(r'real\s+(\d+)m([\d.]+)s', time_str) + if match: + minutes = int(match.group(1)) + seconds = float(match.group(2)) + return minutes * 60 + seconds + return None + +def parse_sim_time(output): + match = re.search(r"Simulation time:\s*([\d.]+)\s*seconds", output) + if match: + return float(match.group(1)) + return None + +def main(): + opticks_home = get_opticks_home() + + parser = argparse.ArgumentParser() + parser.add_argument('-g', '--gdml', type=Path, default=None, help="Path to a custom GDML geometry file (relative to current directory)") + parser.add_argument('-o', '--outpath', type=Path, default=Path('./'), help="Path where the output file will be saved") + args = parser.parse_args() + + # If gdml not provided, use default path relative to OPTICKS_HOME + if args.gdml is None: + gdml_path = opticks_home / 'tests/geom/opticks_raindrop.gdml' + else: + # User-provided path is used as-is (relative to current directory or absolute) + gdml_path = args.gdml + + # run.mac is created in OPTICKS_HOME + run_mac_path = opticks_home / "run.mac" + + geant_file = args.outpath / "timing_geant.txt" + optix_file = args.outpath / "timing_optix.txt" + log_file = args.outpath / "g4logs.txt" + + with open(geant_file, "w") as gfile, open(optix_file, "w") as ofile, open(log_file, "w") as logfile: + for threads in range(50, 0, -1): + times = {} + sim_time_true = None + for flag in ['true', 'false']: + # Write run.mac with current flag + with open(run_mac_path, "w") as rm: + rm.write(run_mac_template.format(threads=threads, flag=flag)) + + # Run with time in bash to capture real/user/sys + cmd = f"time simg4oxmt -g {gdml_path} -m {run_mac_path}" + print(f"Running {threads} threads: {cmd}") + result = subprocess.run( + ["bash", "-c", cmd], + capture_output=True, text=True + ) + stdout = result.stdout + stderr = result.stderr + # Save full output to log file + logfile.write(f"\n{'='*60}\n") + logfile.write(f"Threads: {threads}, StackPhotons: {flag}\n") + logfile.write(f"{'='*60}\n") + logfile.write(f"--- STDOUT ---\n{stdout}\n") + logfile.write(f"--- STDERR ---\n{stderr}\n") + logfile.flush() + + # Save simulation time for true run only + if flag == 'true': + sim_time_true = parse_sim_time(stdout + stderr) + if sim_time_true is not None: + ofile.write(f"{threads} {sim_time_true}\n") + ofile.flush() + + # Extract real time + real_match = re.search(r"real\s+\d+m[\d.]+s", stderr) + if real_match: + real_sec = parse_real_time(real_match.group()) + times[flag] = real_sec + else: + print(f"[!] Could not find 'real' time for threads={threads} flag={flag}") + + # Write the difference to timing_geant.txt (true - false) + if 'true' in times and 'false' in times: + diff = times['true'] - times['false'] + gfile.write(f"{threads} {diff}\n") + gfile.flush() + else: + print(f"[!] Missing times for threads={threads}") + + print("Done.") + +if __name__ == '__main__': + main() diff --git a/pyproject.toml b/pyproject.toml index e3d82cc3aa..7ef5bc28e8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,6 +14,7 @@ dependencies = [ generate-input-photons = "optiphy.tools.generate_input_photons:main" plot-csg = "optiphy.tools.plot_csg:main" serve-path = "optiphy.tools.serve_path:main" +run-performance = "optiphy.tools.run_performance:main" [build-system] requires = ["setuptools", "wheel"] diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index cccbfa8258..2ed4381026 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -32,6 +32,14 @@ target_include_directories(simg4ox PRIVATE $ ) +# simg4ox runs Geant4 and OptiX simulations +add_executable(simg4oxmt simg4oxmt.cpp g4appmt.h) +target_link_libraries(simg4oxmt gphox) +target_include_directories(simg4oxmt PRIVATE + $ + $ +) + # simtox creates a numpy file with initial photons for simulation add_executable(simtox simtox.cpp) @@ -42,7 +50,7 @@ target_include_directories(simtox PRIVATE $ ) -install(TARGETS consgeo simg4ox simtox gphox +install(TARGETS consgeo simg4ox simg4oxmt simtox gphox EXPORT ${PROJECT_NAME}Targets RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} diff --git a/src/g4appmt.h b/src/g4appmt.h new file mode 100644 index 0000000000..b8d3d62199 --- /dev/null +++ b/src/g4appmt.h @@ -0,0 +1,523 @@ +#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 "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; + + EventAction(SEvt *sev) : sev(sev) + { + } + + void BeginOfEventAction(const G4Event *event) override + { + } + + void EndOfEventAction(const G4Event *event) override + { + } +}; + +struct RunAction : G4UserRunAction +{ + RunAction() + { + } + + 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; + } + } +}; + +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); + } + } + + 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++) + { + if ((*procPost)[i3]->GetProcessName() == "Cerenkov") + { + 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; + } + + G4Cerenkov *proc = (G4Cerenkov *)(*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); + + const G4Event *event = G4EventManager::GetEventManager()->GetConstCurrentEvent(); + if (!event) + return; // Always check for null + G4int eventid = event->GetEventID(); + + unsigned int num_hits = 0; + + if (num_hits > 0) + { + G4HCtable *hctable = G4SDManager::GetSDMpointer()->GetHCtable(); + for (G4int i = 0; i < hctable->entries(); ++i) + { + std::string sdn = hctable->GetSDname(i); + std::size_t found = sdn.find("PhotonDetector"); + if (found != std::string::npos) + { + std::cout << "PhotonDetector: " << sdn << std::endl; + PhotonSD *aSD = + (PhotonSD *)G4SDManager::GetSDMpointer()->FindSensitiveDetector(sdn); + } + } + } + } + } + } + } + } +}; + +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()), + 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/src/simg4oxmt.cpp b/src/simg4oxmt.cpp new file mode 100644 index 0000000000..1c2b5b8e98 --- /dev/null +++ b/src/simg4oxmt.cpp @@ -0,0 +1,125 @@ +#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 "g4appmt.h" + +#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_); + } +}; + +int main(int argc, char **argv) +{ + + long seed = static_cast(time(nullptr)); + CLHEP::HepRandom::setTheSeed(seed); + G4cout << "Random seed set to: " << seed << G4endl; + OPTICKS_LOG(argc, argv); + + argparse::ArgumentParser program("simg4ox", "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("whether to open an interactive window with a viewer") + .flag() + .store_into(interactive); + + try + { + program.parse_args(argc, argv); + } + catch (const exception &err) + { + cerr << err.what() << endl; + cerr << program; + exit(EXIT_FAILURE); + } + + // Configure Geant4 + // The physics list must be instantiated before other user actions + G4VModularPhysicsList *physics = new FTFP_BERT; + physics->RegisterPhysics(new G4OpticalPhysics); + + 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/tests/geom/opticks_raindrop.gdml b/tests/geom/opticks_raindrop.gdml index b40f3fc47b..b42151b253 100644 --- a/tests/geom/opticks_raindrop.gdml +++ b/tests/geom/opticks_raindrop.gdml @@ -1,23 +1,16 @@ - + - - - - - - - - - + + + + + @@ -50,8 +43,6 @@ - - @@ -122,8 +113,8 @@ - - + + @@ -131,7 +122,7 @@ - + @@ -151,8 +142,8 @@ - - + + @@ -177,20 +168,14 @@ - - - - + + - - - - - + diff --git a/tests/geom/pfrich_min_FINAL.gdml b/tests/geom/pfrich_min_FINAL.gdml new file mode 100644 index 0000000000..9dcb4fb0f1 --- /dev/null +++ b/tests/geom/pfrich_min_FINAL.gdml @@ -0,0 +1,2677 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +