eic-opticks/dd4hepplugins/OpticsEvent.cc

0001 #include "OpticsEvent.hh"
0002
0003 #include <DD4hep/InstanceCount.h>
0004 #include <DDG4/Factories.h>
0005 #include <DDG4/Geant4Context.h>
0006 #include <DDG4/Geant4Data.h>
0007 #include <DDG4/Geant4HitCollection.h>
0008 #include <DDG4/Geant4SensDetAction.h>
0009
0010 #include <G4Event.hh>
0011
0012 #include <chrono>
0013 #include <cuda_runtime.h>
0014
0015 #include <G4CXOpticks.hh>
0016 #include <NP.hh>
0017 #include <NPFold.h>
0018 #include <SComp.h>
0019 #include <SEvt.hh>
0020 #include <map>
0021 #include <sphoton.h>
0022
0023 namespace ddeicopticks
0024 {
0025 //---------------------------------------------------------------------------//
0026 OpticsEvent::OpticsEvent(dd4hep::sim::Geant4Context *ctxt, std::string const &name)
0027     : dd4hep::sim::Geant4EventAction(ctxt, name)
0028 {
0029     dd4hep::InstanceCount::increment(this);
0030     declareProperty("Verbose", verbose_);
0031     declareProperty("PhotonThreshold", photon_threshold_);
0032 }
0033
0034 //---------------------------------------------------------------------------//
0035 OpticsEvent::~OpticsEvent()
0036 {
0037     dd4hep::InstanceCount::decrement(this);
0038 }
0039
0040 //---------------------------------------------------------------------------//
0041 void OpticsEvent::begin(G4Event const *event)
0042 {
0043     int eventID = event->GetEventID();
0044
0045     if (verbose_ > 0)
0046         info("OpticsEvent::begin -- event #%d", eventID);
0047
0048     // In batch mode, only start a new SEvt at the beginning of each batch.
0049     // Skipping beginOfEvent between events lets gensteps accumulate
0050     // (SEvt::clear_output preserves gensteps by design).
0051     if (!batch_begun_)
0052     {
0053         SEvt::CreateOrReuse_EGPU();
0054         SEvt *sev = SEvt::Get_EGPU();
0055         if (sev)
0056             sev->beginOfEvent(eventID);
0057         batch_begun_ = true;
0058     }
0059 }
0060
0061 //---------------------------------------------------------------------------//
0062 void OpticsEvent::end(G4Event const *event)
0063 {
0064     int eventID = event->GetEventID();
0065
0066     G4CXOpticks *gx = G4CXOpticks::Get();
0067     if (!gx)
0068     {
0069         error("OpticsEvent::end -- G4CXOpticks not initialized");
0070         return;
0071     }
0072
0073     SEvt *sev = SEvt::Get_EGPU();
0074     if (!sev)
0075     {
0076         error("OpticsEvent::end -- no EGPU SEvt instance");
0077         return;
0078     }
0079
0080     int64_t num_genstep = sev->getNumGenstepFromGenstep();
0081     int64_t num_photon = sev->getNumPhotonFromGenstep();
0082
0083     if (verbose_ > 0 || num_genstep > 0)
0084     {
0085         info("Event #%d: %lld gensteps, %lld photons accumulated", eventID, static_cast<long long>(num_genstep),
0086              static_cast<long long>(num_photon));
0087     }
0088
0089     // Batch mode: keep accumulating until threshold reached
0090     if (photon_threshold_ > 0 && num_photon < photon_threshold_)
0091         return;
0092
0093     if (num_genstep > 0)
0094     {
0095         auto sim_t0 = std::chrono::high_resolution_clock::now();
0096         gx->simulate(eventID, /*reset=*/false);
0097         cudaDeviceSynchronize();
0098         auto sim_t1 = std::chrono::high_resolution_clock::now();
0099         double simulate_ms = std::chrono::duration<double, std::milli>(sim_t1 - sim_t0).count();
0100
0101         unsigned num_hit = sev->getNumHit();
0102         info("OPTICKS_GPU_TIME event=%d ms=%.3f photons=%lld hits=%u", eventID, simulate_ms,
0103              static_cast<long long>(num_photon), num_hit);
0104
0105         // Debug: dump photon flag and boundary statistics (verbose >= 2)
0106         if (verbose_ >= 2)
0107         {
0108             NPFold *tf = sev->topfold;
0109             const NP *photon = tf ? tf->get("photon") : nullptr;
0110             if (photon)
0111             {
0112                 int np = photon->shape[0];
0113                 const sphoton *pp = (const sphoton *)photon->cvalues<float>();
0114                 std::map<unsigned, int> flag_counts;
0115                 std::map<unsigned, int> boundary_counts;
0116                 for (int i = 0; i < np; i++)
0117                 {
0118                     flag_counts[pp[i].flagmask]++;
0119                     boundary_counts[pp[i].boundary()]++;
0120                 }
0121                 for (auto &[f, c] : flag_counts)
0122                     info("  flagmask 0x%08x : %d photons", f, c);
0123                 for (auto &[b, c] : boundary_counts)
0124                     info("  boundary %u : %d photons", b, c);
0125             }
0126         }
0127
0128         // Inject hits only in per-event mode; batch mode cannot map
0129         // hits back to individual events.
0130         if (photon_threshold_ == 0 && num_hit > 0)
0131             injectHits(event, sev, num_hit);
0132
0133         sev->endOfEvent(eventID);
0134         gx->reset(eventID);
0135     }
0136     else
0137     {
0138         if (verbose_ > 0)
0139             info("Event #%d: no gensteps, skipping GPU simulation", eventID);
0140         if (photon_threshold_ == 0)
0141             sev->endOfEvent(eventID);
0142     }
0143
0144     batch_begun_ = false;
0145 }
0146
0147 //---------------------------------------------------------------------------//
0148 void OpticsEvent::injectHits(G4Event const *event, SEvt *sev, unsigned num_hit)
0149 {
0150     using dd4hep::sim::Geant4HitCollection;
0151     using dd4hep::sim::Geant4SensDetSequences;
0152     using dd4hep::sim::Geant4Tracker;
0153
0154     int eventID = event->GetEventID();
0155
0156     Geant4SensDetSequences &sens = context()->sensitiveActions();
0157     auto const &seqs = sens.sequences();
0158
0159     if (seqs.empty())
0160     {
0161         warning("Event #%d: no sensitive detectors registered -- "
0162                 "call setupTracker() in steering script",
0163                 eventID);
0164         return;
0165     }
0166
0167     if (seqs.size() > 1)
0168     {
0169         warning("Event #%d: %zu sensitive detectors registered -- "
0170                 "hit routing by sensor identity not yet implemented, "
0171                 "injecting into first SD only",
0172                 eventID, seqs.size());
0173     }
0174
0175     // Inject into the first SD with a valid collection.
0176     // TODO: route hits to the correct SD based on sensor identity
0177     // when multiple sensitive detectors are present.
0178     for (auto const &[det_name, seq] : seqs)
0179     {
0180         Geant4HitCollection *coll = seq->collection(0);
0181         if (!coll)
0182             continue;
0183
0184         for (unsigned i = 0; i < num_hit; i++)
0185         {
0186             sphoton ph;
0187             sev->getHit(ph, i);
0188             coll->add(createTrackerHit(ph));
0189         }
0190
0191         info("Event #%d: injected %u hits into '%s' collection", eventID, num_hit, det_name.c_str());
0192         break;
0193     }
0194 }
0195
0196 //---------------------------------------------------------------------------//
0197 dd4hep::sim::Geant4Tracker::Hit *OpticsEvent::createTrackerHit(sphoton const &ph)
0198 {
0199     using dd4hep::sim::Geant4Tracker;
0200
0201     auto *hit = new Geant4Tracker::Hit();
0202
0203     hit->position = {ph.pos.x, ph.pos.y, ph.pos.z};
0204     hit->momentum = {ph.mom.x, ph.mom.y, ph.mom.z};
0205     hit->length = ph.wavelength;
0206     hit->energyDeposit = 0;
0207     hit->cellID = ph.pmtid();
0208
0209     hit->truth.trackID = ph.index;
0210     hit->truth.pdgID = 0;
0211     hit->truth.deposit = 0;
0212     hit->truth.time = ph.time;
0213     hit->truth.length = ph.wavelength;
0214     hit->truth.x = ph.pos.x;
0215     hit->truth.y = ph.pos.y;
0216     hit->truth.z = ph.pos.z;
0217     hit->truth.px = ph.mom.x;
0218     hit->truth.py = ph.mom.y;
0219     hit->truth.pz = ph.mom.z;
0220
0221     return hit;
0222 }
0223
0224 //---------------------------------------------------------------------------//
0225 } // namespace ddeicopticks
0226
0227 DECLARE_GEANT4ACTION_NS(ddeicopticks, OpticsEvent)