|
| 1 | +#include "OpticsEvent.hh" |
| 2 | + |
| 3 | +#include <DD4hep/InstanceCount.h> |
| 4 | +#include <DDG4/Factories.h> |
| 5 | +#include <DDG4/Geant4Data.h> |
| 6 | +#include <DDG4/Geant4HitCollection.h> |
| 7 | +#include <DDG4/Geant4SensDetAction.h> |
| 8 | +#include <DDG4/Geant4Context.h> |
| 9 | + |
| 10 | +#include <G4Event.hh> |
| 11 | + |
| 12 | +#include <chrono> |
| 13 | + |
| 14 | +#include <G4CXOpticks.hh> |
| 15 | +#include <SEvt.hh> |
| 16 | +#include <SComp.h> |
| 17 | +#include <sphoton.h> |
| 18 | +#include <NP.hh> |
| 19 | +#include <NPFold.h> |
| 20 | +#include <map> |
| 21 | + |
| 22 | +namespace ddeicopticks |
| 23 | +{ |
| 24 | +//---------------------------------------------------------------------------// |
| 25 | +OpticsEvent::OpticsEvent( dd4hep::sim::Geant4Context* ctxt, |
| 26 | + std::string const& name ) |
| 27 | + : dd4hep::sim::Geant4EventAction( ctxt, name ) |
| 28 | +{ |
| 29 | + dd4hep::InstanceCount::increment( this ); |
| 30 | + declareProperty( "Verbose", verbose_ ); |
| 31 | + declareProperty( "PhotonThreshold", photon_threshold_ ); |
| 32 | +} |
| 33 | + |
| 34 | +//---------------------------------------------------------------------------// |
| 35 | +OpticsEvent::~OpticsEvent() |
| 36 | +{ |
| 37 | + dd4hep::InstanceCount::decrement( this ); |
| 38 | +} |
| 39 | + |
| 40 | +//---------------------------------------------------------------------------// |
| 41 | +void OpticsEvent::begin( G4Event const* event ) |
| 42 | +{ |
| 43 | + int eventID = event->GetEventID(); |
| 44 | + |
| 45 | + if( verbose_ > 0 ) |
| 46 | + info( "OpticsEvent::begin -- event #%d", eventID ); |
| 47 | + |
| 48 | + // In batch mode, only start a new SEvt at the beginning of each batch. |
| 49 | + // Skipping beginOfEvent between events lets gensteps accumulate |
| 50 | + // (SEvt::clear_output preserves gensteps by design). |
| 51 | + if( !batch_begun_ ) |
| 52 | + { |
| 53 | + SEvt::CreateOrReuse_EGPU(); |
| 54 | + SEvt* sev = SEvt::Get_EGPU(); |
| 55 | + if( sev ) |
| 56 | + sev->beginOfEvent( eventID ); |
| 57 | + batch_begun_ = true; |
| 58 | + } |
| 59 | +} |
| 60 | + |
| 61 | +//---------------------------------------------------------------------------// |
| 62 | +void OpticsEvent::end( G4Event const* event ) |
| 63 | +{ |
| 64 | + int eventID = event->GetEventID(); |
| 65 | + |
| 66 | + G4CXOpticks* gx = G4CXOpticks::Get(); |
| 67 | + if( !gx ) |
| 68 | + { |
| 69 | + error( "OpticsEvent::end -- G4CXOpticks not initialized" ); |
| 70 | + return; |
| 71 | + } |
| 72 | + |
| 73 | + SEvt* sev = SEvt::Get_EGPU(); |
| 74 | + if( !sev ) |
| 75 | + { |
| 76 | + error( "OpticsEvent::end -- no EGPU SEvt instance" ); |
| 77 | + return; |
| 78 | + } |
| 79 | + |
| 80 | + int64_t num_genstep = sev->getNumGenstepFromGenstep(); |
| 81 | + int64_t num_photon = sev->getNumPhotonFromGenstep(); |
| 82 | + |
| 83 | + if( verbose_ > 0 || num_genstep > 0 ) |
| 84 | + { |
| 85 | + info( "Event #%d: %lld gensteps, %lld photons accumulated", |
| 86 | + eventID, |
| 87 | + static_cast<long long>( num_genstep ), |
| 88 | + static_cast<long long>( num_photon ) ); |
| 89 | + } |
| 90 | + |
| 91 | + // Batch mode: keep accumulating until threshold reached |
| 92 | + if( photon_threshold_ > 0 && num_photon < photon_threshold_ ) |
| 93 | + return; |
| 94 | + |
| 95 | + if( num_genstep > 0 ) |
| 96 | + { |
| 97 | + auto sim_t0 = std::chrono::high_resolution_clock::now(); |
| 98 | + gx->simulate( eventID, /*reset=*/false ); |
| 99 | + auto sim_t1 = std::chrono::high_resolution_clock::now(); |
| 100 | + double simulate_ms = |
| 101 | + std::chrono::duration<double, std::milli>( sim_t1 - sim_t0 ).count(); |
| 102 | + |
| 103 | + unsigned num_hit = sev->getNumHit(); |
| 104 | + info( "OPTICKS_GPU_TIME event=%d ms=%.3f photons=%lld hits=%u", |
| 105 | + eventID, simulate_ms, |
| 106 | + static_cast<long long>( num_photon ), num_hit ); |
| 107 | + |
| 108 | + // Debug: dump photon flag and boundary statistics (verbose >= 2) |
| 109 | + if( verbose_ >= 2 ) |
| 110 | + { |
| 111 | + NPFold* tf = sev->topfold; |
| 112 | + const NP* photon = tf ? tf->get( "photon" ) : nullptr; |
| 113 | + if( photon ) |
| 114 | + { |
| 115 | + int np = photon->shape[0]; |
| 116 | + const sphoton* pp = |
| 117 | + (const sphoton*)photon->cvalues<float>(); |
| 118 | + std::map<unsigned, int> flag_counts; |
| 119 | + std::map<unsigned, int> boundary_counts; |
| 120 | + for( int i = 0; i < np; i++ ) |
| 121 | + { |
| 122 | + flag_counts[pp[i].flagmask]++; |
| 123 | + boundary_counts[pp[i].boundary()]++; |
| 124 | + } |
| 125 | + for( auto& [f, c] : flag_counts ) |
| 126 | + info( " flagmask 0x%08x : %d photons", f, c ); |
| 127 | + for( auto& [b, c] : boundary_counts ) |
| 128 | + info( " boundary %u : %d photons", b, c ); |
| 129 | + } |
| 130 | + } |
| 131 | + |
| 132 | + // Inject hits only in per-event mode; batch mode cannot map |
| 133 | + // hits back to individual events. |
| 134 | + if( photon_threshold_ == 0 && num_hit > 0 ) |
| 135 | + injectHits( event, sev, num_hit ); |
| 136 | + |
| 137 | + sev->endOfEvent( eventID ); |
| 138 | + gx->reset( eventID ); |
| 139 | + } |
| 140 | + else |
| 141 | + { |
| 142 | + if( verbose_ > 0 ) |
| 143 | + info( "Event #%d: no gensteps, skipping GPU simulation", eventID ); |
| 144 | + if( photon_threshold_ == 0 ) |
| 145 | + sev->endOfEvent( eventID ); |
| 146 | + } |
| 147 | + |
| 148 | + batch_begun_ = false; |
| 149 | +} |
| 150 | + |
| 151 | +//---------------------------------------------------------------------------// |
| 152 | +void OpticsEvent::injectHits( G4Event const* event, |
| 153 | + SEvt* sev, |
| 154 | + unsigned num_hit ) |
| 155 | +{ |
| 156 | + using dd4hep::sim::Geant4HitCollection; |
| 157 | + using dd4hep::sim::Geant4SensDetSequences; |
| 158 | + using dd4hep::sim::Geant4Tracker; |
| 159 | + |
| 160 | + int eventID = event->GetEventID(); |
| 161 | + |
| 162 | + Geant4SensDetSequences& sens = context()->sensitiveActions(); |
| 163 | + auto const& seqs = sens.sequences(); |
| 164 | + |
| 165 | + if( seqs.empty() ) |
| 166 | + { |
| 167 | + warning( "Event #%d: no sensitive detectors registered -- " |
| 168 | + "call setupTracker() in steering script", |
| 169 | + eventID ); |
| 170 | + return; |
| 171 | + } |
| 172 | + |
| 173 | + if( seqs.size() > 1 ) |
| 174 | + { |
| 175 | + warning( "Event #%d: %zu sensitive detectors registered -- " |
| 176 | + "hit routing by sensor identity not yet implemented, " |
| 177 | + "injecting into first SD only", |
| 178 | + eventID, seqs.size() ); |
| 179 | + } |
| 180 | + |
| 181 | + // Inject into the first SD with a valid collection. |
| 182 | + // TODO: route hits to the correct SD based on sensor identity |
| 183 | + // when multiple sensitive detectors are present. |
| 184 | + for( auto const& [det_name, seq] : seqs ) |
| 185 | + { |
| 186 | + Geant4HitCollection* coll = seq->collection( 0 ); |
| 187 | + if( !coll ) |
| 188 | + continue; |
| 189 | + |
| 190 | + for( unsigned i = 0; i < num_hit; i++ ) |
| 191 | + { |
| 192 | + sphoton ph; |
| 193 | + sev->getHit( ph, i ); |
| 194 | + coll->add( createTrackerHit( ph ) ); |
| 195 | + } |
| 196 | + |
| 197 | + info( "Event #%d: injected %u hits into '%s' collection", |
| 198 | + eventID, num_hit, det_name.c_str() ); |
| 199 | + break; |
| 200 | + } |
| 201 | +} |
| 202 | + |
| 203 | +//---------------------------------------------------------------------------// |
| 204 | +dd4hep::sim::Geant4Tracker::Hit* |
| 205 | +OpticsEvent::createTrackerHit( sphoton const& ph ) |
| 206 | +{ |
| 207 | + using dd4hep::sim::Geant4Tracker; |
| 208 | + |
| 209 | + auto* hit = new Geant4Tracker::Hit(); |
| 210 | + |
| 211 | + hit->position = { ph.pos.x, ph.pos.y, ph.pos.z }; |
| 212 | + hit->momentum = { ph.mom.x, ph.mom.y, ph.mom.z }; |
| 213 | + hit->length = ph.wavelength; |
| 214 | + hit->energyDeposit = 0; |
| 215 | + hit->cellID = ph.identity; |
| 216 | + |
| 217 | + hit->truth.trackID = ph.index; |
| 218 | + hit->truth.pdgID = 0; |
| 219 | + hit->truth.deposit = 0; |
| 220 | + hit->truth.time = ph.time; |
| 221 | + hit->truth.length = ph.wavelength; |
| 222 | + hit->truth.x = ph.pos.x; |
| 223 | + hit->truth.y = ph.pos.y; |
| 224 | + hit->truth.z = ph.pos.z; |
| 225 | + hit->truth.px = ph.mom.x; |
| 226 | + hit->truth.py = ph.mom.y; |
| 227 | + hit->truth.pz = ph.mom.z; |
| 228 | + |
| 229 | + return hit; |
| 230 | +} |
| 231 | + |
| 232 | +//---------------------------------------------------------------------------// |
| 233 | +} // namespace ddeicopticks |
| 234 | + |
| 235 | +DECLARE_GEANT4ACTION_NS( ddeicopticks, OpticsEvent ) |
0 commit comments