diff --git a/.gitignore b/.gitignore index b25c15b81..5be058f46 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,4 @@ *~ +dd4hepplugins/examples/drich-dev/ +dd4hepplugins/examples/calibrations/ +dd4hepplugins/examples/fieldmaps/ diff --git a/CMakeLists.txt b/CMakeLists.txt index 2917e4003..0825c1059 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -39,6 +39,16 @@ add_subdirectory(g4cx) add_subdirectory(g4cx/tests) add_subdirectory(src) +# Optional DD4hep integration plugins +find_package(DD4hep QUIET COMPONENTS DDG4 DDCore) +if(DD4hep_FOUND) + message(STATUS "DD4hep found -- building dd4hepplugins") + add_subdirectory(dd4hepplugins) + add_subdirectory(dd4hepplugins/examples) +else() + message(STATUS "DD4hep not found -- skipping dd4hepplugins") +endif() + # Export configs include(CMakePackageConfigHelpers) diff --git a/dd4hepplugins/CMakeLists.txt b/dd4hepplugins/CMakeLists.txt new file mode 100644 index 000000000..a3a971b90 --- /dev/null +++ b/dd4hepplugins/CMakeLists.txt @@ -0,0 +1,63 @@ +#------------------------------- -*- cmake -*- -------------------------------# +# eic-opticks DD4hep integration plugins +# +# Builds DD4hep action plugins that integrate eic-opticks GPU-accelerated +# optical photon simulation into DD4hep/Geant4. +# +# Works both as part of the top-level eic-opticks build (in-tree) and as a +# standalone project against an installed eic-opticks. +# +# Usage from a DD4hep steering file: +# OpticsRun -- initializes/finalizes G4CXOpticks per run +# OpticsEvent -- triggers GPU simulation per event +#----------------------------------------------------------------------------# + +# Detect standalone vs in-tree build +if(NOT TARGET G4CX) + cmake_minimum_required(VERSION 3.18 FATAL_ERROR) + project(ddeicopticks VERSION 0.1.0 LANGUAGES CXX) + find_package(DD4hep REQUIRED COMPONENTS DDG4 DDCore) + find_package(eic-opticks REQUIRED) + find_package(Geant4 REQUIRED) + set(_g4cx eic-opticks::G4CX) + set(_u4 eic-opticks::U4) + set(_sysrap eic-opticks::SysRap) +else() + # In-tree: DD4hep already found by parent, targets use local names + find_package(Geant4 REQUIRED) + set(_g4cx G4CX) + set(_u4 U4) + set(_sysrap SysRap) +endif() + +dd4hep_set_compiler_flags() + +set(CMAKE_LIBRARY_OUTPUT_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}") +set(LIBRARY_OUTPUT_PATH "${CMAKE_CURRENT_BINARY_DIR}") + +dd4hep_add_plugin(ddeicopticks + SOURCES + OpticsRun.cc + OpticsEvent.cc + OpticsSteppingAction.cc + USES + DD4hep::DDG4 + DD4hep::DDCore + ${_g4cx} + ${_u4} + ${_sysrap} +) +# STANDALONE changes class export macros in eic-opticks headers; +# only needed when building outside the main eic-opticks tree. +if(NOT TARGET G4CX) + target_compile_definitions(ddeicopticks PRIVATE STANDALONE) +endif() + +install(TARGETS ddeicopticks + LIBRARY DESTINATION lib +) + +install(FILES + ${CMAKE_CURRENT_BINARY_DIR}/${CMAKE_SHARED_MODULE_PREFIX}ddeicopticks.components + DESTINATION lib +) diff --git a/dd4hepplugins/DD4hepSensorIdentifier.hh b/dd4hepplugins/DD4hepSensorIdentifier.hh new file mode 100644 index 000000000..a70d98970 --- /dev/null +++ b/dd4hepplugins/DD4hepSensorIdentifier.hh @@ -0,0 +1,93 @@ +#pragma once +/** +DD4hepSensorIdentifier.hh +=========================== + +Custom sensor identifier for DD4hep geometries. + +Unlike the default U4SensorIdentifierDefault which relies on +GLOBAL_SENSOR_BOUNDARY_LIST env var for non-instanced geometries, +this implementation directly checks G4VSensitiveDetector on volumes. + +This works for DD4hep geometries where sensitive detectors are +explicitly set via DetElement::setSensitiveDetector(). +**/ + +#include +#include + +#include "G4PVPlacement.hh" +#include "G4VSensitiveDetector.hh" + +#include "U4SensorIdentifier.h" + +struct DD4hepSensorIdentifier : public U4SensorIdentifier +{ + int level = 0; + int counter = 0; // auto-increment sensor ID (1-based; 0 means "not a sensor" in opticks) + + void setLevel(int _level) override + { + level = _level; + } + + /** + getGlobalIdentity + ------------------- + Checks if the physical volume has a G4VSensitiveDetector attached. + Returns a unique 1-based sensor_id, or -1 if not sensitive. + + Note: opticks treats sensor_id == 0 as "not a sensor", so IDs must be >= 1. + PV copy numbers are not reliable (e.g. dRICH SiPMs all have copyNo=0). + **/ + int getGlobalIdentity(const G4VPhysicalVolume *pv, const G4VPhysicalVolume * /*ppv*/) override + { + if (!pv) + return -1; + + const G4LogicalVolume *lv = pv->GetLogicalVolume(); + G4VSensitiveDetector *sd = lv->GetSensitiveDetector(); + + if (!sd) + return -1; + + int sensor_id = ++counter; // 1-based unique ID + + if (level > 0) + std::cout << "DD4hepSensorIdentifier::getGlobalIdentity" << " sensor_id " << sensor_id << " sd " + << sd->GetName() << " pv " << pv->GetName() << std::endl; + + return sensor_id; + } + + /** + getInstanceIdentity + --------------------- + Same as default: recursively search for G4VSensitiveDetector + within the instance subtree. + **/ + int getInstanceIdentity(const G4VPhysicalVolume *instance_outer_pv) const override + { + if (!instance_outer_pv) + return -1; + + std::vector sdpv; + FindSD_r(sdpv, instance_outer_pv, 0); + + if (sdpv.empty()) + return -1; + + const G4PVPlacement *pvp = dynamic_cast(instance_outer_pv); + return pvp ? pvp->GetCopyNo() : 0; + } + + static void FindSD_r(std::vector &sdpv, const G4VPhysicalVolume *pv, int depth) + { + const G4LogicalVolume *lv = pv->GetLogicalVolume(); + G4VSensitiveDetector *sd = lv->GetSensitiveDetector(); + if (sd) + sdpv.push_back(pv); + for (size_t i = 0; i < size_t(lv->GetNoDaughters()); i++) + FindSD_r(sdpv, lv->GetDaughter(i), depth + 1); + } +}; diff --git a/dd4hepplugins/OpticsEvent.cc b/dd4hepplugins/OpticsEvent.cc new file mode 100644 index 000000000..5de1fc8e0 --- /dev/null +++ b/dd4hepplugins/OpticsEvent.cc @@ -0,0 +1,227 @@ +#include "OpticsEvent.hh" + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace ddeicopticks +{ +//---------------------------------------------------------------------------// +OpticsEvent::OpticsEvent(dd4hep::sim::Geant4Context *ctxt, std::string const &name) + : dd4hep::sim::Geant4EventAction(ctxt, name) +{ + dd4hep::InstanceCount::increment(this); + declareProperty("Verbose", verbose_); + declareProperty("PhotonThreshold", photon_threshold_); +} + +//---------------------------------------------------------------------------// +OpticsEvent::~OpticsEvent() +{ + dd4hep::InstanceCount::decrement(this); +} + +//---------------------------------------------------------------------------// +void OpticsEvent::begin(G4Event const *event) +{ + int eventID = event->GetEventID(); + + if (verbose_ > 0) + info("OpticsEvent::begin -- event #%d", eventID); + + // In batch mode, only start a new SEvt at the beginning of each batch. + // Skipping beginOfEvent between events lets gensteps accumulate + // (SEvt::clear_output preserves gensteps by design). + if (!batch_begun_) + { + SEvt::CreateOrReuse_EGPU(); + SEvt *sev = SEvt::Get_EGPU(); + if (sev) + sev->beginOfEvent(eventID); + batch_begun_ = true; + } +} + +//---------------------------------------------------------------------------// +void OpticsEvent::end(G4Event const *event) +{ + int eventID = event->GetEventID(); + + G4CXOpticks *gx = G4CXOpticks::Get(); + if (!gx) + { + error("OpticsEvent::end -- G4CXOpticks not initialized"); + return; + } + + SEvt *sev = SEvt::Get_EGPU(); + if (!sev) + { + error("OpticsEvent::end -- no EGPU SEvt instance"); + return; + } + + int64_t num_genstep = sev->getNumGenstepFromGenstep(); + int64_t num_photon = sev->getNumPhotonFromGenstep(); + + if (verbose_ > 0 || num_genstep > 0) + { + info("Event #%d: %lld gensteps, %lld photons accumulated", eventID, static_cast(num_genstep), + static_cast(num_photon)); + } + + // Batch mode: keep accumulating until threshold reached + if (photon_threshold_ > 0 && num_photon < photon_threshold_) + return; + + if (num_genstep > 0) + { + auto sim_t0 = std::chrono::high_resolution_clock::now(); + gx->simulate(eventID, /*reset=*/false); + cudaDeviceSynchronize(); + auto sim_t1 = std::chrono::high_resolution_clock::now(); + double simulate_ms = std::chrono::duration(sim_t1 - sim_t0).count(); + + unsigned num_hit = sev->getNumHit(); + info("OPTICKS_GPU_TIME event=%d ms=%.3f photons=%lld hits=%u", eventID, simulate_ms, + static_cast(num_photon), num_hit); + + // Debug: dump photon flag and boundary statistics (verbose >= 2) + if (verbose_ >= 2) + { + NPFold *tf = sev->topfold; + const NP *photon = tf ? tf->get("photon") : nullptr; + if (photon) + { + int np = photon->shape[0]; + const sphoton *pp = (const sphoton *)photon->cvalues(); + std::map flag_counts; + std::map boundary_counts; + for (int i = 0; i < np; i++) + { + flag_counts[pp[i].flagmask]++; + boundary_counts[pp[i].boundary()]++; + } + for (auto &[f, c] : flag_counts) + info(" flagmask 0x%08x : %d photons", f, c); + for (auto &[b, c] : boundary_counts) + info(" boundary %u : %d photons", b, c); + } + } + + // Inject hits only in per-event mode; batch mode cannot map + // hits back to individual events. + if (photon_threshold_ == 0 && num_hit > 0) + injectHits(event, sev, num_hit); + + sev->endOfEvent(eventID); + gx->reset(eventID); + } + else + { + if (verbose_ > 0) + info("Event #%d: no gensteps, skipping GPU simulation", eventID); + if (photon_threshold_ == 0) + sev->endOfEvent(eventID); + } + + batch_begun_ = false; +} + +//---------------------------------------------------------------------------// +void OpticsEvent::injectHits(G4Event const *event, SEvt *sev, unsigned num_hit) +{ + using dd4hep::sim::Geant4HitCollection; + using dd4hep::sim::Geant4SensDetSequences; + using dd4hep::sim::Geant4Tracker; + + int eventID = event->GetEventID(); + + Geant4SensDetSequences &sens = context()->sensitiveActions(); + auto const &seqs = sens.sequences(); + + if (seqs.empty()) + { + warning("Event #%d: no sensitive detectors registered -- " + "call setupTracker() in steering script", + eventID); + return; + } + + if (seqs.size() > 1) + { + warning("Event #%d: %zu sensitive detectors registered -- " + "hit routing by sensor identity not yet implemented, " + "injecting into first SD only", + eventID, seqs.size()); + } + + // Inject into the first SD with a valid collection. + // TODO: route hits to the correct SD based on sensor identity + // when multiple sensitive detectors are present. + for (auto const &[det_name, seq] : seqs) + { + Geant4HitCollection *coll = seq->collection(0); + if (!coll) + continue; + + for (unsigned i = 0; i < num_hit; i++) + { + sphoton ph; + sev->getHit(ph, i); + coll->add(createTrackerHit(ph)); + } + + info("Event #%d: injected %u hits into '%s' collection", eventID, num_hit, det_name.c_str()); + break; + } +} + +//---------------------------------------------------------------------------// +dd4hep::sim::Geant4Tracker::Hit *OpticsEvent::createTrackerHit(sphoton const &ph) +{ + using dd4hep::sim::Geant4Tracker; + + auto *hit = new Geant4Tracker::Hit(); + + hit->position = {ph.pos.x, ph.pos.y, ph.pos.z}; + hit->momentum = {ph.mom.x, ph.mom.y, ph.mom.z}; + hit->length = ph.wavelength; + hit->energyDeposit = 0; + hit->cellID = ph.pmtid(); + + hit->truth.trackID = ph.index; + hit->truth.pdgID = 0; + hit->truth.deposit = 0; + hit->truth.time = ph.time; + hit->truth.length = ph.wavelength; + hit->truth.x = ph.pos.x; + hit->truth.y = ph.pos.y; + hit->truth.z = ph.pos.z; + hit->truth.px = ph.mom.x; + hit->truth.py = ph.mom.y; + hit->truth.pz = ph.mom.z; + + return hit; +} + +//---------------------------------------------------------------------------// +} // namespace ddeicopticks + +DECLARE_GEANT4ACTION_NS(ddeicopticks, OpticsEvent) diff --git a/dd4hepplugins/OpticsEvent.hh b/dd4hepplugins/OpticsEvent.hh new file mode 100644 index 000000000..b7ca518c1 --- /dev/null +++ b/dd4hepplugins/OpticsEvent.hh @@ -0,0 +1,50 @@ +#pragma once + +#include +#include +#include +#include + +class SEvt; +struct sphoton; + +namespace ddeicopticks +{ +//---------------------------------------------------------------------------// +/*! + * DDG4 action plugin for eic-opticks event-level GPU simulation. + * + * At begin-of-event: prepares GPU event buffer (SEvt). + * At end-of-event: triggers GPU optical photon simulation via + * G4CXOpticks::simulate(), retrieves hits, injects them into DD4hep + * hit collections, and resets for next event. + * + * Requires setupTracker() to be called in the steering script so that + * DD4hep creates hit collections for the sensitive detectors. + * + * Properties: + * - Verbose (default: 0) -- verbosity level + */ +class OpticsEvent final : public dd4hep::sim::Geant4EventAction +{ + public: + OpticsEvent(dd4hep::sim::Geant4Context *ctxt, std::string const &name); + + void begin(G4Event const *event) final; + void end(G4Event const *event) final; + + protected: + DDG4_DEFINE_ACTION_CONSTRUCTORS(OpticsEvent); + ~OpticsEvent() final; + + private: + void injectHits(G4Event const *event, SEvt *sev, unsigned num_hit); + static dd4hep::sim::Geant4Tracker::Hit *createTrackerHit(sphoton const &ph); + + int verbose_{0}; + int64_t photon_threshold_{0}; ///< 0 = simulate per-event, >0 = batch until N photons + bool batch_begun_{false}; +}; + +//---------------------------------------------------------------------------// +} // namespace ddeicopticks diff --git a/dd4hepplugins/OpticsRun.cc b/dd4hepplugins/OpticsRun.cc new file mode 100644 index 000000000..5032e4eee --- /dev/null +++ b/dd4hepplugins/OpticsRun.cc @@ -0,0 +1,155 @@ +#include "OpticsRun.hh" + +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include "DD4hepSensorIdentifier.hh" + +namespace ddeicopticks +{ + +//---------------------------------------------------------------------------// +/*! + * Add Geant4 10.x property name aliases for Opticks compatibility. + * + * Opticks GPU code looks up scintillation spectra by Geant4 10.x names + * (FASTCOMPONENT, SLOWCOMPONENT) while DD4hep + Geant4 11.x uses + * SCINTILLATIONCOMPONENT1/2. This adds the old names pointing to the + * same data so both sides find what they need. + */ +static void addOpticksPropertyAliases() +{ + static const std::pair aliases[] = { + {"SCINTILLATIONCOMPONENT1", "FASTCOMPONENT"}, + {"SCINTILLATIONCOMPONENT2", "SLOWCOMPONENT"}, + }; + + auto const *matTable = G4Material::GetMaterialTable(); + for (auto const *mat : *matTable) + { + auto *mpt = mat->GetMaterialPropertiesTable(); + if (!mpt) + continue; + + for (auto const &[g4_11_name, g4_10_name] : aliases) + { + auto *prop = mpt->GetProperty(g4_11_name); + if (prop && !mpt->GetProperty(g4_10_name)) + { + bool createNewKey = true; + mpt->AddProperty(g4_10_name, prop, createNewKey); + } + } + } +} + +//---------------------------------------------------------------------------// +OpticsRun::OpticsRun(dd4hep::sim::Geant4Context *ctxt, std::string const &name) + : dd4hep::sim::Geant4RunAction(ctxt, name) +{ + dd4hep::InstanceCount::increment(this); + declareProperty("SaveGeometry", save_geometry_); +} + +//---------------------------------------------------------------------------// +OpticsRun::~OpticsRun() +{ + dd4hep::InstanceCount::decrement(this); +} + +//---------------------------------------------------------------------------// +void OpticsRun::begin(G4Run const *run) +{ + G4VPhysicalVolume *world = context()->world(); + if (!world) + { + except("OpticsRun: world volume is null at begin-of-run"); + return; + } + + info("Initializing G4CXOpticks geometry (run #%d)", run->GetRunID()); + + // Add Geant4 10.x scintillation property aliases for Opticks GPU + addOpticksPropertyAliases(); + + SEvt::CreateOrReuse(SEvt::EGPU); + + // Register DD4hep-aware sensor identifier before geometry translation. + // Unlike the default which requires GLOBAL_SENSOR_BOUNDARY_LIST env var, + // this checks G4VSensitiveDetector directly on volumes. + static DD4hepSensorIdentifier dd4hep_sid; + G4CXOpticks::SetSensorIdentifier(&dd4hep_sid); + + bool hasDevice = SEventConfig::HasDevice(); + info("HasDevice=%s, IntegrationMode=%d", hasDevice ? "YES" : "NO", SEventConfig::IntegrationMode()); + G4CXOpticks::SetGeometry(world); + + if (save_geometry_) + { + info("Saving Opticks geometry to disk"); + G4CXOpticks::SaveGeometry(); + } + + // Log boundary count + { + const SSim *sim = SSim::Get(); + if (sim && sim->get_bnd()) + info("Boundary table: %zu boundaries", sim->get_bnd()->names.size()); + } + + info("G4CXOpticks geometry initialized successfully"); +} + +//---------------------------------------------------------------------------// +void OpticsRun::end(G4Run const *run) +{ + // Flush any remaining batched gensteps (from PhotonThreshold mode) + SEvt *sev = SEvt::Get_EGPU(); + if (sev) + { + int64_t num_genstep = sev->getNumGenstepFromGenstep(); + int64_t num_photon = sev->getNumPhotonFromGenstep(); + if (num_genstep > 0) + { + G4CXOpticks *gx = G4CXOpticks::Get(); + if (gx) + { + int eventID = run->GetNumberOfEvent() > 0 ? run->GetNumberOfEvent() - 1 : 0; + info("Flushing %lld remaining photons from %lld gensteps", static_cast(num_photon), + static_cast(num_genstep)); + + auto t0 = std::chrono::high_resolution_clock::now(); + gx->simulate(eventID, /*reset=*/false); + auto t1 = std::chrono::high_resolution_clock::now(); + double ms = std::chrono::duration(t1 - t0).count(); + + unsigned num_hit = sev->getNumHit(); + info("OPTICKS_GPU_TIME event=%d ms=%.3f photons=%lld hits=%u", eventID, ms, + static_cast(num_photon), num_hit); + + sev->endOfEvent(eventID); + gx->reset(eventID); + } + } + } + + info("Finalizing G4CXOpticks (run #%d)", run->GetRunID()); + G4CXOpticks::Finalize(); +} + +//---------------------------------------------------------------------------// +} // namespace ddeicopticks + +DECLARE_GEANT4ACTION_NS(ddeicopticks, OpticsRun) diff --git a/dd4hepplugins/OpticsRun.hh b/dd4hepplugins/OpticsRun.hh new file mode 100644 index 000000000..a9fec1e4b --- /dev/null +++ b/dd4hepplugins/OpticsRun.hh @@ -0,0 +1,38 @@ +#pragma once + +#include +#include +#include + +namespace ddeicopticks +{ +//---------------------------------------------------------------------------// +/*! + * DDG4 action plugin for eic-opticks run-level lifecycle. + * + * At begin-of-run: initializes G4CXOpticks geometry from the Geant4 world + * volume, translating the geometry to OptiX acceleration structures on GPU. + * + * At end-of-run: finalizes G4CXOpticks, releasing GPU resources. + * + * Properties: + * - SaveGeometry (default: false) -- save Opticks geometry to disk + */ +class OpticsRun final : public dd4hep::sim::Geant4RunAction +{ + public: + OpticsRun(dd4hep::sim::Geant4Context *ctxt, std::string const &name); + + void begin(G4Run const *run) final; + void end(G4Run const *run) final; + + protected: + DDG4_DEFINE_ACTION_CONSTRUCTORS(OpticsRun); + ~OpticsRun() final; + + private: + bool save_geometry_{false}; +}; + +//---------------------------------------------------------------------------// +} // namespace ddeicopticks diff --git a/dd4hepplugins/OpticsSteppingAction.cc b/dd4hepplugins/OpticsSteppingAction.cc new file mode 100644 index 000000000..655242b17 --- /dev/null +++ b/dd4hepplugins/OpticsSteppingAction.cc @@ -0,0 +1,144 @@ +#include "OpticsSteppingAction.hh" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace +{ +G4Mutex genstep_mutex = G4MUTEX_INITIALIZER; +} + +namespace ddeicopticks +{ +//---------------------------------------------------------------------------// +OpticsSteppingAction::OpticsSteppingAction(dd4hep::sim::Geant4Context *ctxt, std::string const &name) + : dd4hep::sim::Geant4SteppingAction(ctxt, name) +{ + dd4hep::InstanceCount::increment(this); + declareProperty("Verbose", verbose_); +} + +//---------------------------------------------------------------------------// +OpticsSteppingAction::~OpticsSteppingAction() +{ + dd4hep::InstanceCount::decrement(this); +} + +//---------------------------------------------------------------------------// +void OpticsSteppingAction::operator()(const G4Step *step, G4SteppingManager *mgr) +{ + G4SteppingManager *fpSteppingManager = mgr; + + const G4Track *track = step->GetTrack(); + + G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus(); + if (stepStatus == fAtRestDoItProc) + return; + + // Skip optical photons — they don't produce gensteps + if (track->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) + return; + + G4ProcessVector *procPost = fpSteppingManager->GetfPostStepDoItVector(); + size_t nproc = fpSteppingManager->GetMAXofPostStepLoops(); + + if (verbose_ > 1) + { + std::string procs; + for (size_t i = 0; i < nproc; i++) + if ((*procPost)[i]) + procs += (*procPost)[i]->GetProcessName() + " "; + info(" nproc=%zu processes=[%s]", nproc, procs.c_str()); + } + + for (size_t i = 0; i < nproc; i++) + { + if (!(*procPost)[i]) + continue; + + const G4String &procName = (*procPost)[i]->GetProcessName(); + + if (procName == "Cerenkov") + { + const G4DynamicParticle *aParticle = track->GetDynamicParticle(); + G4double charge = aParticle->GetDefinition()->GetPDGCharge(); + const G4Material *mat = track->GetMaterial(); + G4MaterialPropertiesTable *mpt = mat->GetMaterialPropertiesTable(); + if (!mpt) + continue; + + G4MaterialPropertyVector *Rindex = mpt->GetProperty(kRINDEX); + if (!Rindex || Rindex->GetVectorLength() == 0) + continue; + + G4Cerenkov *cer = static_cast((*procPost)[i]); + G4int numPhotons = cer->GetNumPhotons(); + if (numPhotons <= 0) + continue; + + G4double Pmin = Rindex->Energy(0); + G4double Pmax = Rindex->GetMaxEnergy(); + G4double nMax = Rindex->GetMaxValue(); + G4double beta1 = step->GetPreStepPoint()->GetBeta(); + G4double beta2 = step->GetPostStepPoint()->GetBeta(); + G4double beta = (beta1 + beta2) * 0.5; + G4double BetaInverse = 1.0 / beta; + G4double maxCos = BetaInverse / nMax; + G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos); + G4double MeanNumberOfPhotons1 = cer->GetAverageNumberOfPhotons(charge, beta1, mat, Rindex); + G4double MeanNumberOfPhotons2 = cer->GetAverageNumberOfPhotons(charge, beta2, mat, Rindex); + + { + G4AutoLock lock(&genstep_mutex); + U4::CollectGenstep_G4Cerenkov_modified(track, step, numPhotons, BetaInverse, Pmin, Pmax, maxCos, + maxSin2, MeanNumberOfPhotons1, MeanNumberOfPhotons2); + } + + if (verbose_ > 0) + info("Cerenkov genstep: %d photons", numPhotons); + } + else if (procName == "Scintillation") + { + G4Scintillation *scint = static_cast((*procPost)[i]); + G4int numPhotons = scint->GetNumPhotons(); + if (numPhotons <= 0) + continue; + + const G4Material *mat = track->GetMaterial(); + G4MaterialPropertiesTable *mpt = mat->GetMaterialPropertiesTable(); + if (!mpt) + continue; + + G4double scintTime = 0.0; + if (mpt->ConstPropertyExists(kSCINTILLATIONTIMECONSTANT1)) + scintTime = mpt->GetConstProperty(kSCINTILLATIONTIMECONSTANT1); + + { + G4AutoLock lock(&genstep_mutex); + U4::CollectGenstep_DsG4Scintillation_r4695(track, step, numPhotons, + /*scnt=*/1, scintTime); + } + + if (verbose_ > 0) + info("Scintillation genstep: %d photons", numPhotons); + } + } +} + +//---------------------------------------------------------------------------// +} // namespace ddeicopticks + +DECLARE_GEANT4ACTION_NS(ddeicopticks, OpticsSteppingAction) diff --git a/dd4hepplugins/OpticsSteppingAction.hh b/dd4hepplugins/OpticsSteppingAction.hh new file mode 100644 index 000000000..0439b0f90 --- /dev/null +++ b/dd4hepplugins/OpticsSteppingAction.hh @@ -0,0 +1,34 @@ +#pragma once + +#include +#include +#include + +namespace ddeicopticks +{ +//---------------------------------------------------------------------------// +/*! + * DDG4 stepping action that intercepts standard Geant4 Cerenkov and + * Scintillation processes and collects gensteps for GPU simulation. + * + * This follows the same approach as eic-opticks GPUCerenkov example: + * use standard G4Cerenkov / G4Scintillation, then read back the photon + * count via GetNumPhotons() and pack gensteps for the GPU. + */ +class OpticsSteppingAction final : public dd4hep::sim::Geant4SteppingAction +{ + public: + OpticsSteppingAction(dd4hep::sim::Geant4Context *ctxt, std::string const &name); + + void operator()(const G4Step *step, G4SteppingManager *mgr) final; + + protected: + DDG4_DEFINE_ACTION_CONSTRUCTORS(OpticsSteppingAction); + ~OpticsSteppingAction() final; + + private: + int verbose_{0}; +}; + +//---------------------------------------------------------------------------// +} // namespace ddeicopticks diff --git a/dd4hepplugins/examples/CMakeLists.txt b/dd4hepplugins/examples/CMakeLists.txt new file mode 100644 index 000000000..326573783 --- /dev/null +++ b/dd4hepplugins/examples/CMakeLists.txt @@ -0,0 +1,3 @@ +find_package(DD4hep REQUIRED COMPONENTS DDCore) +dd4hep_set_compiler_flags() +dd4hep_add_plugin(RaindropGeo SOURCES geometry/Raindrop_geo.cpp USES DD4hep::DDCore) diff --git a/dd4hepplugins/examples/benchmark_drich.py b/dd4hepplugins/examples/benchmark_drich.py new file mode 100644 index 000000000..3b50364ed --- /dev/null +++ b/dd4hepplugins/examples/benchmark_drich.py @@ -0,0 +1,393 @@ +#!/usr/bin/env python3 +""" +Benchmark: GPU speedup for optical photon simulation on the EPIC dRICH detector. + +Runs three configurations with the full dRICH geometry (100 pi+ at 10 GeV): + + cpu -- Optical photons generated AND tracked on CPU by Geant4 + baseline -- Optical photons generated but NOT tracked (SetStackPhotons=false) + gpu -- Optical photons simulated on GPU via eic-opticks + +Speedup = CPU_optical_time / GPU_simulate_time + where CPU_optical_time = T(cpu) - T(baseline) + and GPU_simulate_time = wall time of G4CXOpticks::simulate() + +Prerequisites: + - Spack environment activated with DD4hep, Geant4, eic-opticks + - `spack load epic` (sets DETECTOR_PATH to EPIC geometry) + - DD4HEP_LIBRARY_PATH includes /opt/local/lib (eic-opticks plugins) + +Usage: + python3 benchmark_drich.py # all modes, 10 events + python3 benchmark_drich.py --events 100 --photon-threshold 5000000 + python3 benchmark_drich.py --mode gpu --events 10 +""" +import argparse +import json +import math +import os +import re +import subprocess +import sys +import time + + +_SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) + + +def find_compact_file(geometry="1sector"): + """Locate the dRICH compact XML. + + geometry: '1sector' (default) uses the bundled single-sector XML, + 'full' uses EPIC's epic_drich_only.xml (all 6 sectors). + """ + if geometry == "1sector": + compact = os.path.join(_SCRIPT_DIR, "geometry", "epic_drich_1sector.xml") + os.environ["DRICH_1SECTOR_XML"] = os.path.join( + _SCRIPT_DIR, "geometry", "drich_1sector.xml") + else: + detector_path = os.environ.get("DETECTOR_PATH", "") + if not detector_path: + print("ERROR: DETECTOR_PATH not set. Run: spack load epic", + file=sys.stderr) + sys.exit(1) + compact = os.path.join(detector_path, "epic_drich_only.xml") + + if not os.path.exists(compact): + print(f"ERROR: {compact} not found", file=sys.stderr) + sys.exit(1) + + return compact + + +def run_single_mode(mode, num_events, photon_threshold=0, multiplicity=100, + geometry="1sector"): + """Run one benchmark mode inside the current process.""" + import cppyy + import DDG4 + from g4units import GeV + + cppyy.include("G4OpticalParameters.hh") + from cppyy.gbl import G4OpticalParameters + + compact = find_compact_file(geometry) + + # ---- kernel & geometry ------------------------------------------------ + kernel = DDG4.Kernel() + print(f"Loading geometry: {compact}") + kernel.loadGeometry(str("file:" + compact)) + + geant4 = DDG4.Geant4(kernel) + geant4.printDetectors() + geant4.setupUI(typ="tcsh", vis=False, ui=False) + + # ---- particle gun: pi+ at 10 GeV toward dRICH sector 0 --------------- + # dRICH acceptance: eta 1.6-3.5 => theta 3.5-22.8 deg + # ideal theta: eta=2.0 => ~15.4 deg (full rings in one sector) + # 1-sector geometry: sector 0 sensors are at phi ~ -167 deg + # Full geometry: sector 0 at same phi; other sectors fill the ring + eta = 2.0 + theta = 2.0 * math.atan(math.exp(-eta)) # ~15.4 deg + phi = math.radians(-167.0) if geometry == "1sector" else 0.0 + geant4.setupGun( + "Gun", + particle="pi+", + energy=10 * GeV, + position=(0, 0, 0), + isotrop=False, + direction=( + math.sin(theta) * math.cos(phi), + math.sin(theta) * math.sin(phi), + math.cos(theta), + ), + multiplicity=multiplicity, + ) + + # ---- physics ---------------------------------------------------------- + geant4.setupPhysics("QGSP_BERT") + ph = DDG4.PhysicsList(kernel, "Geant4PhysicsList/OpticalPhys") + ph.addPhysicsConstructor(str("G4OpticalPhysics")) + kernel.physicsList().adopt(ph) + + # ---- mode-specific setup ---------------------------------------------- + if mode == "cpu": + # Track optical photons on CPU with proper optical tracker + seq, act = geant4.setupDetector("DRICH", "Geant4OpticalTrackerAction") + filt = DDG4.Filter( + kernel, "ParticleSelectFilter/OpticalPhotonSelector" + ) + filt.particle = "opticalphoton" + seq.adopt(filt) + + elif mode == "baseline": + # Need a tracker so DD4hep is happy, but no photons will arrive + geant4.setupTracker("DRICH") + + elif mode == "gpu": + # Opticks env vars for hit gathering + os.environ.setdefault("OPTICKS_EVENT_MODE", "Minimal") + os.environ.setdefault("OPTICKS_INTEGRATION_MODE", "1") + + # Block CPU G4Step hits; only GPU-injected hits pass + seq, act = geant4.setupTracker("DRICH") + filt = DDG4.Filter(kernel, "EnergyDepositMinimumCut/Block") + filt.Cut = 1e12 + seq.adopt(filt) + + # eic-opticks DDG4 plugins + stepping = DDG4.SteppingAction( + kernel, "OpticsSteppingAction/OpticsStep1" + ) + stepping.Verbose = 0 + kernel.steppingAction().adopt(stepping) + + run_action = DDG4.RunAction(kernel, "OpticsRun/OpticsRun1") + run_action.SaveGeometry = False + kernel.runAction().adopt(run_action) + + evt_action = DDG4.EventAction(kernel, "OpticsEvent/OpticsEvt1") + evt_action.Verbose = 1 + if photon_threshold > 0: + evt_action.PhotonThreshold = photon_threshold + kernel.eventAction().adopt(evt_action) + + # ---- configure & initialize ------------------------------------------- + kernel.NumEvents = num_events + kernel.configure() + + # Disable photon stacking BEFORE initialize so G4Cerenkov reads it + if mode in ("baseline", "gpu"): + G4OpticalParameters.Instance().SetCerenkovStackPhotons(False) + G4OpticalParameters.Instance().SetScintStackPhotons(False) + + kernel.initialize() + + # BoundaryInvokeSD AFTER initialize (runtime parameter) + if mode == "cpu": + G4OpticalParameters.Instance().SetBoundaryInvokeSD(True) + + # ---- run & time ------------------------------------------------------- + t0 = time.perf_counter() + kernel.run() + t1 = time.perf_counter() + wall_s = t1 - t0 + + kernel.terminate() + + result = { + "mode": mode, + "events": num_events, + "wall_s": round(wall_s, 4), + "per_event_ms": round(wall_s / num_events * 1000, 2), + "photon_threshold": photon_threshold, + "multiplicity": multiplicity, + "geometry": geometry, + } + print(f"BENCHMARK_RESULT:{json.dumps(result)}", flush=True) + return result + + +# --------------------------------------------------------------------------- +# "all" mode — runs each config as a subprocess, then computes speedup +# --------------------------------------------------------------------------- + +def _run_subprocess(mode, num_events, photon_threshold=0, multiplicity=100, + geometry="1sector"): + """Run a single mode in a child process, return (result_dict, raw_output).""" + script = os.path.abspath(__file__) + cmd = [ + sys.executable, script, + "--mode", mode, + "--events", str(num_events), + "--photon-threshold", str(photon_threshold), + "--multiplicity", str(multiplicity), + "--geometry", geometry, + ] + proc = subprocess.run( + cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True + ) + output = proc.stdout + + if proc.returncode != 0: + print(f"ERROR: subprocess exited with code {proc.returncode}", + file=sys.stderr) + return None, output + + result = None + for line in output.splitlines(): + if line.startswith("BENCHMARK_RESULT:"): + result = json.loads(line[len("BENCHMARK_RESULT:"):]) + break + + return result, output + + +def _parse_gpu_times(output): + """Extract per-event GPU simulate() times from C++ info output.""" + times = [] + for line in output.splitlines(): + m = re.search(r"OPTICKS_GPU_TIME event=(\d+) ms=([\d.]+)", line) + if m: + times.append((int(m.group(1)), float(m.group(2)))) + return times + + +def _parse_gpu_photons(output): + """Extract total photon count from GPU output.""" + total = 0 + for line in output.splitlines(): + m = re.search(r"OPTICKS_GPU_TIME.*photons=(\d+)", line) + if m: + total += int(m.group(1)) + return total + + +def run_all(num_events, photon_threshold=0, multiplicity=100, geometry="1sector"): + """Run all three modes and print speedup analysis.""" + results = {} + gpu_event_times = [] + gpu_total_photons = 0 + + for mode in ("baseline", "cpu", "gpu"): + pt = photon_threshold if mode == "gpu" else 0 + label = { + "baseline": "baseline (no photon tracking)", + "cpu": "cpu (photons tracked on CPU)", + "gpu": "gpu (photons on GPU via eic-opticks)", + }[mode] + extra = f", threshold={pt}" if pt > 0 else "" + print(f"\n{'='*60}") + print(f" {label} -- {num_events} events{extra}") + phi_str = "phi=-167 (sector 0)" if geometry == "1sector" else "phi=0" + print(f" {multiplicity} pi+ at 10 GeV, eta=2.0 (theta~15.4deg), {phi_str}") + print(f"{'='*60}") + + result, output = _run_subprocess(mode, num_events, pt, multiplicity, + geometry) + + if result is None: + print(f"ERROR: mode={mode} failed. Output:\n{output[-3000:]}") + sys.exit(1) + + results[mode] = result + + # Show selected output lines + for line in output.splitlines(): + if line.startswith("BENCHMARK_RESULT:"): + continue + if any(kw in line for kw in ( + "OPTICKS_GPU_TIME", "gensteps", "hits from GPU", + "Detected photons", + )): + print(f" {line.strip()}") + + print(f" Wall time: {result['wall_s']:.3f} s " + f"({result['per_event_ms']:.1f} ms/event)") + + if mode == "gpu": + gpu_event_times = _parse_gpu_times(output) + gpu_total_photons = _parse_gpu_photons(output) + + # ------------------------------------------------------------------ + # Summary + # ------------------------------------------------------------------ + T_b = results["baseline"]["wall_s"] + T_c = results["cpu"]["wall_s"] + T_g = results["gpu"]["wall_s"] + cpu_optical = T_c - T_b + gpu_overhead = T_g - T_b + + photons_per_event = gpu_total_photons / num_events if num_events > 0 else 0 + + print(f"\n{'='*60}") + print(f" RESULTS ({num_events} events, {multiplicity} pi+ @ 10 GeV in dRICH)") + print(f"{'='*60}") + print(f"\n {'Mode':<12} {'Total (s)':<12} {'Per event (ms)'}") + print(f" {'-'*40}") + for m in ("baseline", "cpu", "gpu"): + r = results[m] + print(f" {m:<12} {r['wall_s']:<12.3f} {r['per_event_ms']:.1f}") + + print(f"\n CPU optical photon tracking: {cpu_optical:.3f} s" + f" ({cpu_optical / num_events * 1000:.1f} ms/event)") + print(f" GPU total overhead: {gpu_overhead:.3f} s" + f" ({gpu_overhead / num_events * 1000:.1f} ms/event)") + + if gpu_total_photons > 0: + print(f"\n Total Cerenkov photons: {gpu_total_photons:,}" + f" (~{photons_per_event:,.0f}/event)") + + if gpu_event_times: + all_ms = [t for _, t in gpu_event_times] + total_ms = sum(all_ms) + print(f"\n GPU simulate() calls: {len(all_ms)}") + # Show first 5 and last 2 if many events + show = gpu_event_times + if len(show) > 10: + show = gpu_event_times[:5] + [("...", None)] + gpu_event_times[-2:] + for item in show: + if item[1] is None: + print(f" ...") + else: + evt_id, ms = item + tag = " (includes OptiX warmup)" if evt_id == 0 else "" + print(f" event {evt_id:>3d}: {ms:>8.1f} ms{tag}") + print(f" {'total':>9s}: {total_ms:>8.1f} ms") + + if len(all_ms) > 1: + warm_ms = all_ms[1:] + avg_warm = sum(warm_ms) / len(warm_ms) + print(f" avg (excl. first): {avg_warm:>5.1f} ms") + + if total_ms > 0: + speedup = cpu_optical / (total_ms / 1000) + print(f"\n >>> SPEEDUP (CPU optical / GPU simulate): {speedup:.1f}x <<<") + + if len(all_ms) > 1 and avg_warm > 0: + speedup_warm = (cpu_optical / num_events * 1000) / avg_warm + print(f" >>> SPEEDUP (per-event, excl. warmup): {speedup_warm:.1f}x <<<") + + if gpu_overhead > 0: + speedup_total = cpu_optical / gpu_overhead + print(f" >>> SPEEDUP (CPU optical / GPU overhead): {speedup_total:.1f}x <<<") + + print(f"{'='*60}\n") + + +# --------------------------------------------------------------------------- +if __name__ == "__main__": + parser = argparse.ArgumentParser( + description="Benchmark GPU vs CPU optical photon speedup on EPIC dRICH" + ) + parser.add_argument( + "--mode", + choices=["cpu", "baseline", "gpu", "all"], + default="all", + help="cpu=photons on CPU, baseline=no photon tracking, " + "gpu=photons on GPU, all=run all three", + ) + parser.add_argument( + "--events", type=int, default=10, + help="number of events (default: 10)", + ) + parser.add_argument( + "--photon-threshold", type=int, default=0, + help="GPU batch mode: accumulate photons across events and simulate " + "when this count is reached (default: 0 = per-event)", + ) + parser.add_argument( + "--multiplicity", type=int, default=100, + help="number of pi+ per event (default: 100)", + ) + parser.add_argument( + "--geometry", choices=["1sector", "full"], default="1sector", + help="1sector=single dRICH sector (default), full=all 6 sectors", + ) + args = parser.parse_args() + + if args.mode == "all": + run_all(args.events, args.photon_threshold, args.multiplicity, + args.geometry) + else: + run_single_mode(args.mode, args.events, args.photon_threshold, + args.multiplicity, args.geometry) diff --git a/dd4hepplugins/examples/benchmark_gpu_speedup.py b/dd4hepplugins/examples/benchmark_gpu_speedup.py new file mode 100644 index 000000000..15dadc2a6 --- /dev/null +++ b/dd4hepplugins/examples/benchmark_gpu_speedup.py @@ -0,0 +1,306 @@ +#!/usr/bin/env python3 +""" +Benchmark: quantify GPU speedup for optical photon simulation. + +Runs three configurations of the raindrop geometry (10 MeV e- in water): + + cpu -- Optical photons generated AND tracked on CPU by Geant4 + baseline -- Optical photons generated but NOT tracked (SetStackPhotons=false) + gpu -- Optical photons simulated on GPU via eic-opticks + +Speedup = CPU_optical_time / GPU_simulate_time + where CPU_optical_time = T(cpu) - T(baseline) + and GPU_simulate_time = wall time of G4CXOpticks::simulate() + (includes genstep upload, OptiX launch, + cudaDeviceSynchronize, hit download) + +Usage: + python3 benchmark_gpu_speedup.py # all modes, 10 events + python3 benchmark_gpu_speedup.py --events 20 # more events + python3 benchmark_gpu_speedup.py --mode cpu # single mode + python3 benchmark_gpu_speedup.py --mode gpu --events 5 # single mode +""" +import argparse +import json +import os +import re +import subprocess +import sys +import time + + +_SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) + + +def run_single_mode(mode, num_events, photon_threshold=0): + """Run one benchmark mode inside the current process.""" + import cppyy + import DDG4 + from g4units import MeV + + cppyy.include("G4OpticalParameters.hh") + from cppyy.gbl import G4OpticalParameters + + # ---- kernel & geometry ------------------------------------------------ + kernel = DDG4.Kernel() + compact = os.path.join(_SCRIPT_DIR, "geometry", "raindrop_dd4hep.xml") + if not os.path.exists(compact): + print(f"ERROR: {compact} not found", file=sys.stderr) + sys.exit(1) + kernel.loadGeometry(str("file:" + compact)) + + geant4 = DDG4.Geant4(kernel) + geant4.setupUI(typ="tcsh", vis=False, ui=False) + geant4.setupGun( + "Gun", + particle="e-", + energy=10 * MeV, + position=(0, 0, 0), + isotrop=False, + direction=(1.0, 0.0, 0.0), + multiplicity=1, + ) + + # ---- physics ---------------------------------------------------------- + geant4.setupPhysics("QGSP_BERT") + ph = DDG4.PhysicsList(kernel, "Geant4PhysicsList/OpticalPhys") + ph.addPhysicsConstructor(str("G4OpticalPhysics")) + kernel.physicsList().adopt(ph) + + # ---- mode-specific setup ---------------------------------------------- + if mode == "cpu": + # Track optical photons on CPU + seq, act = geant4.setupDetector( + "Raindrop", "Geant4OpticalTrackerAction" + ) + filt = DDG4.Filter( + kernel, "ParticleSelectFilter/OpticalPhotonSelector" + ) + filt.particle = "opticalphoton" + seq.adopt(filt) + + elif mode == "baseline": + # Need a tracker so DD4hep is happy, but no photons will arrive + geant4.setupTracker("Raindrop") + + elif mode == "gpu": + # Tracker with impossibly high cut -> blocks CPU-side G4Step hits; + # only GPU-injected hits pass. + seq, act = geant4.setupTracker("Raindrop") + filt = DDG4.Filter(kernel, "EnergyDepositMinimumCut/Block") + filt.Cut = 1e12 + seq.adopt(filt) + + # eic-opticks DDG4 plugins + stepping = DDG4.SteppingAction( + kernel, "OpticsSteppingAction/OpticsStep1" + ) + stepping.Verbose = 0 + kernel.steppingAction().adopt(stepping) + + run_action = DDG4.RunAction(kernel, "OpticsRun/OpticsRun1") + run_action.SaveGeometry = False + kernel.runAction().adopt(run_action) + + evt_action = DDG4.EventAction(kernel, "OpticsEvent/OpticsEvt1") + evt_action.Verbose = 1 + if photon_threshold > 0: + evt_action.PhotonThreshold = photon_threshold + kernel.eventAction().adopt(evt_action) + + # ---- configure & initialize ------------------------------------------- + kernel.NumEvents = num_events + kernel.configure() + + # Disable photon stacking BEFORE initialize so G4Cerenkov reads it + # during BuildPhysicsTable. Cerenkov still fires (gensteps collected) + # but secondary photon tracks are not pushed onto the Geant4 stack. + if mode in ("baseline", "gpu"): + G4OpticalParameters.Instance().SetCerenkovStackPhotons(False) + G4OpticalParameters.Instance().SetScintStackPhotons(False) + + kernel.initialize() + + # BoundaryInvokeSD must be set AFTER initialize (runtime parameter). + # Needed so photons detected at a boundary surface actually invoke the SD. + if mode == "cpu": + G4OpticalParameters.Instance().SetBoundaryInvokeSD(True) + + # ---- run & time ------------------------------------------------------- + t0 = time.perf_counter() + kernel.run() + t1 = time.perf_counter() + wall_s = t1 - t0 + + kernel.terminate() + + result = { + "mode": mode, + "events": num_events, + "wall_s": round(wall_s, 4), + "per_event_ms": round(wall_s / num_events * 1000, 2), + "photon_threshold": photon_threshold, + } + print(f"BENCHMARK_RESULT:{json.dumps(result)}", flush=True) + return result + + +# --------------------------------------------------------------------------- +# "all" mode — runs each config as a subprocess, then computes speedup +# --------------------------------------------------------------------------- + +def _run_subprocess(mode, num_events, photon_threshold=0): + """Run a single mode in a child process, return (result_dict, raw_output).""" + script = os.path.abspath(__file__) + cmd = [sys.executable, script, "--mode", mode, "--events", str(num_events), + "--photon-threshold", str(photon_threshold)] + proc = subprocess.run( + cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True + ) + output = proc.stdout + + if proc.returncode != 0: + print(f"ERROR: subprocess exited with code {proc.returncode}", + file=sys.stderr) + return None, output + + # Parse BENCHMARK_RESULT JSON + result = None + for line in output.splitlines(): + if line.startswith("BENCHMARK_RESULT:"): + result = json.loads(line[len("BENCHMARK_RESULT:"):]) + break + + return result, output + + +def _parse_gpu_times(output): + """Extract per-event GPU simulate() times from C++ info output.""" + times = [] + for line in output.splitlines(): + m = re.search(r"OPTICKS_GPU_TIME event=(\d+) ms=([\d.]+)", line) + if m: + times.append((int(m.group(1)), float(m.group(2)))) + return times + + +def run_all(num_events, photon_threshold=0): + """Run all three modes and print speedup analysis.""" + results = {} + gpu_event_times = [] + + for mode in ("baseline", "cpu", "gpu"): + pt = photon_threshold if mode == "gpu" else 0 + label = { + "baseline": "baseline (no photon tracking)", + "cpu": "cpu (photons tracked on CPU)", + "gpu": "gpu (photons on GPU via eic-opticks)", + }[mode] + extra = f", threshold={pt}" if pt > 0 else "" + print(f"\n{'='*60}") + print(f" {label} -- {num_events} events{extra}") + print(f"{'='*60}") + + result, output = _run_subprocess(mode, num_events, pt) + + if result is None: + print(f"ERROR: mode={mode} failed. Output:\n{output[-2000:]}") + sys.exit(1) + + results[mode] = result + + # Show selected output lines + for line in output.splitlines(): + if line.startswith("BENCHMARK_RESULT:"): + continue + # Show event summary lines and GPU timing + if any(kw in line for kw in ( + "OPTICKS_GPU_TIME", "gensteps", "hits from GPU", + "Detected photons", + )): + print(f" {line.strip()}") + + print(f" Wall time: {result['wall_s']:.3f} s " + f"({result['per_event_ms']:.1f} ms/event)") + + if mode == "gpu": + gpu_event_times = _parse_gpu_times(output) + + # ------------------------------------------------------------------ + # Summary + # ------------------------------------------------------------------ + T_b = results["baseline"]["wall_s"] + T_c = results["cpu"]["wall_s"] + T_g = results["gpu"]["wall_s"] + cpu_optical = T_c - T_b + gpu_overhead = T_g - T_b + + print(f"\n{'='*60}") + print(f" RESULTS ({num_events} events, 10 MeV e- in raindrop)") + print(f"{'='*60}") + print(f"\n {'Mode':<12} {'Total (s)':<12} {'Per event (ms)'}") + print(f" {'-'*40}") + for m in ("baseline", "cpu", "gpu"): + r = results[m] + print(f" {m:<12} {r['wall_s']:<12.3f} {r['per_event_ms']:.1f}") + + print(f"\n CPU optical photon tracking: {cpu_optical:.3f} s" + f" ({cpu_optical / num_events * 1000:.1f} ms/event)") + print(f" GPU total overhead: {gpu_overhead:.3f} s" + f" ({gpu_overhead / num_events * 1000:.1f} ms/event)") + + if gpu_event_times: + all_ms = [t for _, t in gpu_event_times] + total_ms = sum(all_ms) + print(f"\n GPU simulate() per event:") + for evt_id, ms in gpu_event_times: + tag = " (includes OptiX warmup)" if evt_id == 0 else "" + print(f" event {evt_id:>3d}: {ms:>8.1f} ms{tag}") + print(f" {'total':>9s}: {total_ms:>8.1f} ms") + + if len(all_ms) > 1: + warm_ms = all_ms[1:] + avg_warm = sum(warm_ms) / len(warm_ms) + print(f" avg (excl. first): {avg_warm:>5.1f} ms") + + if total_ms > 0: + speedup = cpu_optical / (total_ms / 1000) + print(f"\n >>> SPEEDUP (CPU optical / GPU simulate): {speedup:.1f}x <<<") + + if len(all_ms) > 1 and avg_warm > 0: + speedup_warm = (cpu_optical / num_events * 1000) / avg_warm + print(f" >>> SPEEDUP (per-event, excl. warmup): {speedup_warm:.1f}x <<<") + + if gpu_overhead > 0: + speedup_total = cpu_optical / gpu_overhead + print(f" >>> SPEEDUP (CPU optical / GPU overhead): {speedup_total:.1f}x <<<") + + print(f"{'='*60}\n") + + +# --------------------------------------------------------------------------- +if __name__ == "__main__": + parser = argparse.ArgumentParser( + description="Benchmark GPU vs CPU optical photon speedup" + ) + parser.add_argument( + "--mode", + choices=["cpu", "baseline", "gpu", "all"], + default="all", + help="cpu=photons on CPU, baseline=no photon tracking, " + "gpu=photons on GPU, all=run all three", + ) + parser.add_argument( + "--events", type=int, default=10, help="number of events (default: 10)" + ) + parser.add_argument( + "--photon-threshold", type=int, default=0, + help="GPU batch mode: accumulate photons across events and simulate " + "when this count is reached (default: 0 = per-event)", + ) + args = parser.parse_args() + + if args.mode == "all": + run_all(args.events, args.photon_threshold) + else: + run_single_mode(args.mode, args.events, args.photon_threshold) diff --git a/dd4hepplugins/examples/geometry/Raindrop_geo.cpp b/dd4hepplugins/examples/geometry/Raindrop_geo.cpp new file mode 100644 index 000000000..5ba79aaa7 --- /dev/null +++ b/dd4hepplugins/examples/geometry/Raindrop_geo.cpp @@ -0,0 +1,80 @@ +//========================================================================== +// Raindrop detector constructor for DD4hep +// Based on DD4hep OpNovice example +// +// Geometry: Vacuum(240) > Lead(220) > Air(200) > Water(100) mm +//========================================================================== +#include +#include +#include +#include + +using namespace dd4hep; + +static Ref_t create_raindrop(Detector &description, xml_h e, SensitiveDetector sens) +{ + xml_det_t x_det = e; + std::string name = x_det.nameStr(); + DetElement sdet(name, x_det.id()); + + // Read child elements from XML + xml_det_t x_container = x_det.child(_Unicode(container)); + xml_det_t x_medium = x_det.child(_Unicode(medium)); + xml_det_t x_drop = x_det.child(_Unicode(drop)); + + // Build volumes: nested boxes + Box container_box(x_container.x(), x_container.y(), x_container.z()); + Volume container_vol("CtPMT", container_box, description.material(x_container.attr(_U(material)))); + + Box medium_box(x_medium.x(), x_medium.y(), x_medium.z()); + Volume medium_vol("Md", medium_box, description.material(x_medium.attr(_U(material)))); + + Box drop_box(x_drop.x(), x_drop.y(), x_drop.z()); + Volume drop_vol("Dr", drop_box, description.material(x_drop.attr(_U(material)))); + + // Visualization + if (x_container.hasAttr(_U(vis))) + container_vol.setVisAttributes(description, x_container.visStr()); + if (x_medium.hasAttr(_U(vis))) + medium_vol.setVisAttributes(description, x_medium.visStr()); + if (x_drop.hasAttr(_U(vis))) + drop_vol.setVisAttributes(description, x_drop.visStr()); + + // Mark container (OpLead) as sensitive + container_vol.setSensitiveDetector(sens); + + // Place: drop inside medium, medium inside container + PlacedVolume drop_pv = medium_vol.placeVolume(drop_vol); + drop_pv.addPhysVolID("drop", 1); + + PlacedVolume medium_pv = container_vol.placeVolume(medium_vol); + medium_pv.addPhysVolID("medium", 1); + + // Place container into world + PlacedVolume container_pv = description.pickMotherVolume(sdet).placeVolume(container_vol); + container_pv.addPhysVolID("system", x_det.id()); + sdet.setPlacement(container_pv); + + // Border surface between medium (air) and container (lead) with EFFICIENCY=1.0 + OpticalSurfaceManager surfMgr = description.surfaceManager(); + OpticalSurface surf = surfMgr.opticalSurface("/world/" + name + "#MediumContainerSurf"); + if (surf.isValid()) + { + BorderSurface bsurf = BorderSurface(description, sdet, "MediumContainer", surf, medium_pv, container_pv); + bsurf.isValid(); + printout(INFO, "Raindrop", "Attached border surface MediumContainer (air/lead)"); + } + + // Skin surface on container for Opticks sensor detection + OpticalSurface skinSurf = surfMgr.opticalSurface("/world/" + name + "#ContainerSkinSurf"); + if (skinSurf.isValid()) + { + SkinSurface ssf = SkinSurface(description, sdet, "ContainerSkin", skinSurf, container_vol); + ssf.isValid(); + printout(INFO, "Raindrop", "Attached skin surface ContainerSkin with EFFICIENCY"); + } + + return sdet; +} + +DECLARE_DETELEMENT(DD4hep_Raindrop, create_raindrop) diff --git a/dd4hepplugins/examples/geometry/drich_1sector.xml b/dd4hepplugins/examples/geometry/drich_1sector.xml new file mode 100644 index 000000000..ca9d3b50f --- /dev/null +++ b/dd4hepplugins/examples/geometry/drich_1sector.xml @@ -0,0 +1,362 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +- `DRICH_debug_optics`: 1 = all components become vacuum, except for mirrors; test opticalphotons from IP + 2 = all components become vacuum, except for mirrors and `gasvol`, test charged particles from IP + 3 = all components become vacuum, except for mirrors and sensors, test reflected trajectories' hits + 0 = off +- `DRICH_debug_sector`: 1 = only include one sector (will be set automatically when `DRICH_debug_optics>0`) +- `DRICH_debug_mirror`: 1 = draw full mirror shape for single sector; 0 = off +- `DRICH_debug_sensors`: 1 = draw full sensor sphere for a single sector; 0 = off + + + + + + + + + + + + + +### dRICH: ***d***ual ***R***ing ***I***maging ***Ch***erenkov detector + + + + + + +#### Vessel +- the dRICH vessel is composed of two parts: + - tank: cylindrical region containing most of the detector components + - snout: conical region at the front of the vessel, containing the aerogel +- dimensions: + - `zmin`: z-position of vessel front plane + - `length`: overall z-length of the full vessel + - `snout_length`: length of cone-shaped snout region, housing aerogel + - `rmin0` and `rmin1`: bore radius at front plane and back plane, respectively + - `rmax0` and `rmax1`: outer radius of snout at front plane and snout-back (tank-front) plane, respectively + - `rmax2`: outer radius of tank, the main cylindrical vessel volume + - `nsectors`: number of azimuthal sectors + - `wall_thickness`: thickness of radial walls + - `window_thickness`: thickness of entrance and exit disks + + + + + + +#### Radiator +- radiator is defined in a wedge of azimuthal space, composed of aerogel and a + filter; the filter is applied to the back of the aerogel, so that it separates + the aerogel and gas radiators; an airgap is defined between the aerogel and filter +- dimensions: + - `frontplane`: front of the aerogel, w.r.t. front plane of the vessel envelope + - `rmin` and `rmax`: inner and outer radius (at the front plane; radial bounds are conical) + - `thickness`: radiator thickness, defined separately for aerogel and filter + - `pitch`: controls the angle of the radiator (0=vertical) + + + + + + + + + + +#### Spherical mirror +- spherical mirrors are built from spherical patches, and positioned near the + vessel back plane, separately for each sector +- dimensions: + - `backplane`: the position of the maximum z-plane intersected by the sphere, + w.r.t. the back plane of vessel envelope + - `rmin` and `rmax`: polar angle boundaries + - `phiw`: azimuthal width of one sector + - `thickness` is the radial thickness of the mirror; note that `backplane` is given for the + reflective mirror surface, the inner radius of the sphere + - `focus_tune*` are tuning parameters for the focal region: + - `focus_tune_z` and `focus_tune_x` will move the focal region, with respect + to the sensor sphere center (i.e., set both to zero for focus at the sensor sphere center + (ignoring spherical aberrations effects)) + + + + + +#### Sensors + + + + + +##### Sensor modules: photosensitive surface (pss) + resin base +- SiPM: Silicon Photomultiplier + - PSS: photosensitive surface + - resin: resin substrate (holds the PSS) + - dimensions: + - `side`: side length + - `thickness`: thickness +- PDU: Photodetector Unit (dimensions defined as constants above) + - note: the PDU pitch will determine how many PDUs there are, since the + sensor placement algorithm will try to place as many as it can in the specified + spherical patch below + - Front Services: cooling, heat exchange, etc. + - so far just a box with `side` and `thickness` dimensions + - includes a unique `name`, and `material` and `vis` + - if more than one is defined, they will be layered, in order + - Board: a circuit board (front end board, readout board, etc.) + - dimensions: + - `width`: board side length, paralell to PSS + - `length`: board side length, perpendicular to PSS + - `thickness: board thickness + - `offset`: positioning, with respect to PDU central axis + - Back Services: sockets, cooling, etc. + - so far just a box with `side` and `thickness` dimensions + - includes a unique `name`, and `material` and `vis` + - if more than one is defined, they will be layered, in order + +###### Diagrams +``` + sensor assembly diagram:, where '0' denotes the origin: + + axes: z + +-+--------0--------+-+ | + | | pss | | 0--x + | +-----------------+ | + | resin | + +---------------------+ + +``` +``` +photodetector unit (PDU) assembly diagram: matrix of SiPMs with services + + Top view: 2x2 matrix of SiPMs (2 PDU units shown side-by-side) + ============================= + + -: :- PDU gap size + : : + +-----------+ +-----------+ + | +--+ +--+ | | +--+ +--+ | + | | | | | | | | | | | | + | +--+ +--+ | | +--+ +--+ | + | +--+ +--+ | | +--+ +--+ | + | | | | | | | | | | | | + | +--+ +--+ | | +--+ +--+ | + +-----------+ +-----------+ + : : : : + : : -: :- sensor gap size + : : + -: :- sensor gap size (same) + + Side view: + ========== + + +--------------------+ + | SiPM Matrix | + +--------------------+ + | Cooling, heat | frontservices + | exchange, etc. | + +--------------------+ + || || || || || + || || || || || front-end and + || || || || || readout boards + || || || || || + || || || || || + || + || + +----------------+ + | Sockets, etc. | backservices + +----------------+ + +``` + + + + + + + + + + + + + + + + + + + + + + +##### Sensor sphere +- sensors will be placed on a sphere, using a "disco ball" tiling algorithm; each + sector has its own sensor sphere + - sphere dimensions: + - `centerx` and `centerz`: sphere center, defined w.r.t. vessel front plane, + for the sector on +x axis + - `radius`: radius of the sensor sphere +- sensors will be limited to a patch of the sphere + - patch dimensions: + - `phiw`: defines half the angle between the azimuthal boundaries + - `rmin` and `rmax`: radial cut boundaries + - `zmin`: z-plane cut + + + + + + + + + + + +#### Readout +- segmentation: square matrix of pixels + - `grid_size_x,y`: size of each sensor pixel + - `offset_x,y`: specified such that the `x` and `y` values are unsigned + (except for very rare hits on the sensor sides) + + + + + system:8,sector:3,pdu:12,sipm:6,x:32:-16,y:-16 + + + diff --git a/dd4hepplugins/examples/geometry/epic_drich_1sector.xml b/dd4hepplugins/examples/geometry/epic_drich_1sector.xml new file mode 100644 index 000000000..c5b85b7e3 --- /dev/null +++ b/dd4hepplugins/examples/geometry/epic_drich_1sector.xml @@ -0,0 +1,76 @@ + + + + + + + + + + + + + + EPIC dRICH single sector + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dd4hepplugins/examples/geometry/raindrop_dd4hep.xml b/dd4hepplugins/examples/geometry/raindrop_dd4hep.xml new file mode 100644 index 000000000..ff964eb3f --- /dev/null +++ b/dd4hepplugins/examples/geometry/raindrop_dd4hep.xml @@ -0,0 +1,128 @@ + + + + Minimal optical geometry for eic-opticks DD4hep integration testing. + Matches eic-opticks tests/geom/opticks_raindrop_with_scintillation.gdml + Nested boxes: Vacuum(240) > Lead(220) > Air(200) > Water(100) mm + Water has RINDEX and scintillation. Air has RINDEX. + Border surface between air medium and lead container (MediumContainerSurf). + Skin surface on lead container (ContainerSkinSurf). + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:8,medium:4,drop:4 + + + + diff --git a/dd4hepplugins/examples/test_raindrop_dd4hep_cpu.py b/dd4hepplugins/examples/test_raindrop_dd4hep_cpu.py new file mode 100644 index 000000000..b1cdee366 --- /dev/null +++ b/dd4hepplugins/examples/test_raindrop_dd4hep_cpu.py @@ -0,0 +1,134 @@ +#!/usr/bin/env python3 +""" +CPU-only raindrop optical photon test via DD4hep. + +Uses standard G4OpticalPhysics (no eic-opticks GPU plugins). +Optical photons are tracked on CPU by Geant4 and collected as +Geant4Tracker::Hit via the sensitive detector. + +Compare hit counts with test_raindrop_dd4hep_gpu.py to validate. + +Prerequisites: + - Spack environment activated (ROOT, DD4hep on PYTHONPATH/LD_LIBRARY_PATH) + - DD4hepINSTALL set (for elements.xml lookup) + - libRaindropGeo.so on DD4HEP_LIBRARY_PATH +""" +import os +import sys + +import cppyy +import DDG4 +from g4units import MeV + +cppyy.include("G4OpticalParameters.hh") + +_SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) + + +def run(): + kernel = DDG4.Kernel() + + compact_file = os.path.join(_SCRIPT_DIR, "geometry", "raindrop_dd4hep.xml") + if not os.path.exists(compact_file): + print(f"ERROR: Compact file not found: {compact_file}") + sys.exit(1) + + if "DD4hepINSTALL" not in os.environ: + print("ERROR: DD4hepINSTALL not set. Activate the spack environment first.") + sys.exit(1) + + print(f"Loading geometry: {compact_file}") + kernel.loadGeometry(str("file:" + compact_file)) + + geant4 = DDG4.Geant4(kernel) + geant4.printDetectors() + geant4.setupUI(typ='tcsh', vis=False, ui=False) + + # Same particle gun as the GPU test: e- at 10 MeV inside water box + geant4.setupGun( + "Gun", + particle='e-', + energy=10 * MeV, + position=(0, 0, 0), + isotrop=False, + direction=(1.0, 0.0, 0.0), + multiplicity=1, + ) + + # Register optical photon sensitive detector on lead container. + # Use Geant4OpticalTrackerAction (not default Geant4SimpleTrackerAction + # which has a 1 keV min energy deposit cut -- optical photons are ~eV). + seq, act = geant4.setupDetector('Raindrop', 'Geant4OpticalTrackerAction') + filt = DDG4.Filter(kernel, 'ParticleSelectFilter/OpticalPhotonSelector') + filt.particle = 'opticalphoton' + seq.adopt(filt) + + # Physics: QGSP_BERT + standard Geant4 optical physics (CPU tracking) + geant4.setupPhysics('QGSP_BERT') + ph = DDG4.PhysicsList(kernel, 'Geant4PhysicsList/OpticalPhys') + ph.addPhysicsConstructor(str('G4OpticalPhysics')) + kernel.physicsList().adopt(ph) + + # Save hits to ROOT file + output_file = "/tmp/raindrop_hits_cpu.root" + geant4.setupROOTOutput('RootOutput', output_file) + + # Run 2 events (same as GPU test) + kernel.NumEvents = 2 + kernel.configure() + kernel.initialize() + + # Enable boundary SD invocation so photons detected at the air/lead + # surface (EFFICIENCY=1.0) trigger hits in the sensitive lead volume. + # Geant4 11 defaults BoundaryInvokeSD to false. + from cppyy.gbl import G4OpticalParameters + G4OpticalParameters.Instance().SetBoundaryInvokeSD(True) + + kernel.run() + kernel.terminate() + + # Verify hits + verify_hits(output_file) + + +def verify_hits(root_file): + """Read back detected photon hits and print summary.""" + import ROOT + f = ROOT.TFile.Open(root_file) + if not f or f.IsZombie(): + print(f"ERROR: Cannot open {root_file}") + return + + tree = f.Get("EVENT") + if not tree: + print("ERROR: No EVENT tree in ROOT file") + f.Close() + return + + print(f"\n{'='*60}") + print(f" CPU hit verification from {root_file}") + print(f"{'='*60}") + print(f" Events in file: {tree.GetEntries()}") + + for evt_idx in range(tree.GetEntries()): + tree.GetEntry(evt_idx) + hits = tree.RaindropHits + nhits = hits.size() + + print(f"\n Event {evt_idx}:") + print(f" Detected photons: {nhits}") + for i in range(min(3, nhits)): + h = hits[i] + pos = h.position + mom = h.momentum + print(f" [{i}] pos=({pos.x():.1f}, {pos.y():.1f}, {pos.z():.1f}) " + f"mom=({mom.x():.3f}, {mom.y():.3f}, {mom.z():.3f}) " + f"time={h.truth.time:.2f}ns") + + f.Close() + print(f"\n CPU test complete") + print(f"{'='*60}\n") + + +if __name__ == "__main__": + run() diff --git a/dd4hepplugins/examples/test_raindrop_dd4hep_gpu.py b/dd4hepplugins/examples/test_raindrop_dd4hep_gpu.py new file mode 100644 index 000000000..ad6957295 --- /dev/null +++ b/dd4hepplugins/examples/test_raindrop_dd4hep_gpu.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python3 +""" +Test eic-opticks DD4hep plugins with the raindrop geometry. + +Geometry: Vacuum world > Lead(220mm) > Air(200mm) > Water(100mm) + - Water has RINDEX=1.333 and scintillation properties + - Air has RINDEX=1.0 + - One border surface between water drop and air medium + - Very short volume names (Ct, Md, Dr) -- avoids eic-opticks buffer overflow + +This uses the same geometry as eic-opticks's own raindrop test, +but expressed as DD4hep compact XML instead of GDML. + +Approach: standard G4OpticalPhysics (G4Cerenkov + G4Scintillation) runs on CPU. +OpticsSteppingAction intercepts the processes and collects gensteps. +OpticsEvent triggers GPU simulation and injects hits into DD4hep collections. + +Prerequisites: + - Spack environment activated (ROOT, DD4hep, eic-opticks on PYTHONPATH/LD_LIBRARY_PATH) + - DD4hepINSTALL set (for elements.xml lookup) + - libddeicopticks.so and libRaindropGeo.so on DD4HEP_LIBRARY_PATH +""" +import os +import sys + +import DDG4 +from g4units import GeV, MeV, mm + +_SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) + +def run(): + kernel = DDG4.Kernel() + + compact_file = os.path.join(_SCRIPT_DIR, "geometry", "raindrop_dd4hep.xml") + if not os.path.exists(compact_file): + print(f"ERROR: Compact file not found: {compact_file}") + sys.exit(1) + + if "DD4hepINSTALL" not in os.environ: + print("ERROR: DD4hepINSTALL not set. Activate the spack environment first.") + sys.exit(1) + + print(f"Loading geometry: {compact_file}") + kernel.loadGeometry(str("file:" + compact_file)) + + geant4 = DDG4.Geant4(kernel) + geant4.printDetectors() + + # Batch mode + geant4.setupUI(typ='tcsh', vis=False, ui=False) + + # Particle gun: e- at 10 MeV starting inside water box (50mm half-size) + # Low energy so electron stops inside the water volume + gun = geant4.setupGun( + "Gun", + particle='e-', + energy=10 * MeV, + position=(0, 0, 0), + isotrop=False, + direction=(1.0, 0.0, 0.0), + multiplicity=1, + ) + + # Register sensitive detector -- creates DD4hep hit collection for Raindrop. + # Add EnergyDepositMinimumCut filter to block G4Step hits from charged + # particles; only GPU optical photon hits injected by OpticsEvent pass. + seq, act = geant4.setupTracker('Raindrop') + filt = DDG4.Filter(kernel, 'EnergyDepositMinimumCut/OpticalOnly') + filt.Cut = 1e12 # 1 TeV -- effectively blocks all G4Step hits + seq.adopt(filt) + + # Physics: QGSP_BERT + standard G4OpticalPhysics (G4Cerenkov + G4Scintillation) + geant4.setupPhysics('QGSP_BERT') + ph = DDG4.PhysicsList(kernel, 'Geant4PhysicsList/OpticalPhys') + ph.addPhysicsConstructor(str('G4OpticalPhysics')) + kernel.physicsList().adopt(ph) + + # --- eic-opticks GPU plugins --- + + # SteppingAction: intercepts G4Cerenkov/G4Scintillation and collects gensteps + stepping = DDG4.SteppingAction(kernel, 'OpticsSteppingAction/OpticsStep1') + stepping.Verbose = 1 + kernel.steppingAction().adopt(stepping) + + # RunAction: initializes/finalizes G4CXOpticks geometry + run_action = DDG4.RunAction(kernel, 'OpticsRun/OpticsRun1') + run_action.SaveGeometry = False + kernel.runAction().adopt(run_action) + + # EventAction: triggers GPU simulation, injects hits + evt_action = DDG4.EventAction(kernel, 'OpticsEvent/OpticsEvt1') + evt_action.Verbose = 1 + kernel.eventAction().adopt(evt_action) + + # Save hits to ROOT file via DD4hep's standard output mechanism + output_file = "/tmp/raindrop_hits.root" + geant4.setupROOTOutput('RootOutput', output_file) + + # Run 2 events + kernel.NumEvents = 2 + kernel.configure() + kernel.initialize() + kernel.run() + kernel.terminate() + + # Verify: read back hits from ROOT file + verify_hits(output_file) + +def verify_hits(root_file): + """Read back hits from the ROOT file and print summary.""" + import ROOT + f = ROOT.TFile.Open(root_file) + if not f or f.IsZombie(): + print(f"ERROR: Cannot open {root_file}") + return + + tree = f.Get("EVENT") + if not tree: + print("ERROR: No EVENT tree in ROOT file") + f.Close() + return + + print(f"\n{'='*60}") + print(f" Hit verification from {root_file}") + print(f"{'='*60}") + print(f" Events in file: {tree.GetEntries()}") + + # List branches + branches = [b.GetName() for b in tree.GetListOfBranches()] + print(f" Branches: {branches}") + + hit_branches = [b for b in branches if 'hit' in b.lower() or 'Hit' in b] + if not hit_branches: + hit_branches = [b for b in branches if b not in ('RunNumber', 'EventNumber')] + + for evt_idx in range(tree.GetEntries()): + tree.GetEntry(evt_idx) + print(f"\n Event {evt_idx}:") + for bname in hit_branches: + branch = tree.GetBranch(bname) + leaf = branch.GetLeaf(bname) + if hasattr(tree, bname): + hits = getattr(tree, bname) + nhits = hits.size() if hasattr(hits, 'size') else 0 + print(f" {bname}: {nhits} hits") + # Print first 3 hits + for i in range(min(3, nhits)): + h = hits[i] + pos = h.position + mom = h.momentum + print(f" [{i}] pos=({pos.x():.1f}, {pos.y():.1f}, {pos.z():.1f}) " + f"mom=({mom.x():.3f}, {mom.y():.3f}, {mom.z():.3f}) " + f"wavelength={h.length:.1f}nm " + f"time={h.truth.time:.2f}ns") + + f.Close() + print(f"\n PASS: Hits successfully stored and retrieved via DD4hep") + print(f"{'='*60}\n") + +if __name__ == "__main__": + run() diff --git a/dd4hepplugins/examples/validate_drich.py b/dd4hepplugins/examples/validate_drich.py new file mode 100644 index 000000000..ca33d5edc --- /dev/null +++ b/dd4hepplugins/examples/validate_drich.py @@ -0,0 +1,141 @@ +#!/usr/bin/env python3 +"""Run GPU and CPU with simplified dRICH geometry (no non-optical PDU components).""" +import math, os, sys, numpy as np + +_SCRIPT_DIR = "/tmp/eic-opticks/dd4hepplugins/examples" +_GEOM_DIR = os.path.join(_SCRIPT_DIR, "geometry") +sys.path.insert(0, _SCRIPT_DIR) + +NUM_EVENTS = int(os.environ.get("NUM_EVENTS", "3")) +MULTIPLICITY = int(os.environ.get("MULTIPLICITY", "20")) + + +def setup_kernel(mode): + import cppyy, DDG4 + from g4units import GeV + + cppyy.include("G4OpticalParameters.hh") + from cppyy.gbl import G4OpticalParameters + + compact = os.path.join(_GEOM_DIR, "epic_drich_simple.xml") + if "DRICH_SIMPLE_XML" not in os.environ: + os.environ["DRICH_SIMPLE_XML"] = os.path.join(_GEOM_DIR, "drich_1sector_simple.xml") + + if mode == "gpu": + os.environ["OPTICKS_EVENT_MODE"] = "Minimal" + os.environ["OPTICKS_INTEGRATION_MODE"] = "1" + # Ensure DDG4 can find the plugin + import ctypes + ctypes.CDLL("libddeicopticks.so", ctypes.RTLD_GLOBAL) + else: + os.environ.pop("OPTICKS_EVENT_MODE", None) + os.environ.pop("OPTICKS_INTEGRATION_MODE", None) + + kernel = DDG4.Kernel() + kernel.loadGeometry(str("file:" + compact)) + geant4 = DDG4.Geant4(kernel) + geant4.setupUI(typ="tcsh", vis=False, ui=False) + + eta = 2.0 + theta = 2.0 * math.atan(math.exp(-eta)) + phi = math.radians(-167.0) + geant4.setupGun("Gun", particle="pi+", energy=10*GeV, + position=(0, 0, 0), isotrop=False, + direction=(math.sin(theta)*math.cos(phi), + math.sin(theta)*math.sin(phi), + math.cos(theta)), + multiplicity=MULTIPLICITY) + + geant4.setupPhysics("QGSP_BERT") + ph = DDG4.PhysicsList(kernel, "Geant4PhysicsList/OpticalPhys") + ph.addPhysicsConstructor(str("G4OpticalPhysics")) + kernel.physicsList().adopt(ph) + + # Detector setup needed for both modes (GPU injects hits into DD4hep) + seq, act = geant4.setupDetector("DRICH", "Geant4OpticalTrackerAction") + filt = DDG4.Filter(kernel, "ParticleSelectFilter/OpticalPhotonSelector") + filt.particle = "opticalphoton" + seq.adopt(filt) + + if mode == "gpu": + step_action = DDG4.SteppingAction(kernel, "OpticsSteppingAction/OpticsStep1") + kernel.steppingAction().adopt(step_action) + run_action = DDG4.RunAction(kernel, "OpticsRun/OpticsRun1") + kernel.runAction().adopt(run_action) + # OpticsEvent must be registered BEFORE ROOT output so hits are + # injected into the collection before serialization. + evt_action = DDG4.EventAction(kernel, "OpticsEvent/OpticsEvt1") + kernel.eventAction().adopt(evt_action) + + if mode == "gpu": + geant4.setupROOTOutput("RootOutput", "/tmp/drich_simple_gpu") + else: + geant4.setupROOTOutput("RootOutput", "/tmp/drich_simple_cpu") + + kernel.NumEvents = NUM_EVENTS + kernel.configure() + + # In GPU mode, tell G4Cerenkov/G4Scintillation to NOT push optical + # photon secondaries onto the Geant4 stack. The genstep is still + # computed (numPhotons, BetaInverse, etc.) and collected by + # OpticsSteppingAction, but no CPU photon tracks are created. + if mode == "gpu": + G4OpticalParameters.Instance().SetCerenkovStackPhotons(False) + G4OpticalParameters.Instance().SetScintStackPhotons(False) + + kernel.initialize() + G4OpticalParameters.Instance().SetBoundaryInvokeSD(True) + return kernel + + +def run_gpu(): + kernel = setup_kernel("gpu") + kernel.run() + kernel.terminate() + + # Read hits from npy files + base = os.path.join(os.environ.get("TMP", "/tmp/opticks_out"), + "GEOM", os.environ.get("GEOM", "simple"), + "python3.13", "ALL0_no_opticks_event_name") + total = 0 + if os.path.isdir(base): + for d in sorted(os.listdir(base)): + hp = os.path.join(base, d, "hit.npy") + if os.path.exists(hp): + n = len(np.load(hp)) + total += n + print(f" {d}: {n} hits") + print(f"GPU total: {total} hits") + return total + + +def run_cpu(): + kernel = setup_kernel("cpu") + kernel.run() + kernel.terminate() + + # Extract hits from ROOT + import ROOT + ROOT.gSystem.Load("libDDG4Plugins") + ROOT.gSystem.Load("libDDG4") + f = ROOT.TFile.Open("/tmp/drich_simple_cpu.root") + tree = f.Get("EVENT") + total = 0 + for i in range(tree.GetEntries()): + tree.GetEntry(i) + nhits = tree.DRICHHits.size() + total += nhits + print(f" event {i}: {nhits} hits") + f.Close() + print(f"CPU total: {total} hits") + return total + + +if __name__ == "__main__": + mode = sys.argv[1] if len(sys.argv) > 1 else "gpu" + if mode == "gpu": + run_gpu() + elif mode == "cpu": + run_cpu() + else: + print(f"Usage: {sys.argv[0]} [gpu|cpu]")