eic-opticks/src/GPUPhotonFileSource.h

0001 #include <filesystem>
0002 #include <fstream>
0003 #include <sstream>
0004 #include <string>
0005 #include <vector>
0006
0007 #include "G4Event.hh"
0008 #include "G4GDMLParser.hh"
0009 #include "G4OpticalPhoton.hh"
0010 #include "G4PhysicalConstants.hh"
0011 #include "G4PrimaryParticle.hh"
0012 #include "G4PrimaryVertex.hh"
0013 #include "G4SystemOfUnits.hh"
0014 #include "G4ThreeVector.hh"
0015 #include "G4UserEventAction.hh"
0016 #include "G4UserRunAction.hh"
0017 #include "G4VPhysicalVolume.hh"
0018 #include "G4VUserDetectorConstruction.hh"
0019 #include "G4VUserPrimaryGeneratorAction.hh"
0020
0021 #include "g4cx/G4CXOpticks.hh"
0022 #include "sysrap/NP.hh"
0023 #include "sysrap/SEvt.hh"
0024 #include "sysrap/sphoton.h"
0025
0026 struct DetectorConstruction : G4VUserDetectorConstruction
0027 {
0028     DetectorConstruction(std::filesystem::path gdml_file) : gdml_file_(gdml_file)
0029     {
0030     }
0031
0032     G4VPhysicalVolume *Construct() override
0033     {
0034         parser_.Read(gdml_file_.string(), false);
0035         G4VPhysicalVolume *world = parser_.GetWorldVolume();
0036
0037         G4CXOpticks::SetGeometry(world);
0038
0039         return world;
0040     }
0041
0042   private:
0043     std::filesystem::path gdml_file_;
0044     G4GDMLParser parser_;
0045 };
0046
0047 /// Read photons from a text file.
0048 /// Each line is one photon with 11 space-separated values:
0049 ///   pos_x pos_y pos_z time mom_x mom_y mom_z pol_x pol_y pol_z wavelength
0050 /// Lines starting with '#' are comments and are skipped.
0051 inline std::vector<sphoton> load_photons_txt(const std::filesystem::path &path)
0052 {
0053     std::vector<sphoton> result;
0054     std::ifstream ifs(path);
0055     if (!ifs.is_open())
0056     {
0057         G4cerr << "ERROR: cannot open photon file: " << path << G4endl;
0058         return result;
0059     }
0060
0061     std::string line;
0062     int lineno = 0;
0063     while (std::getline(ifs, line))
0064     {
0065         lineno++;
0066         // skip blank lines and comments
0067         if (line.empty() || line[0] == '#')
0068             continue;
0069
0070         std::istringstream ss(line);
0071         float px, py, pz, t, mx, my, mz, polx, poly, polz, wl;
0072         if (!(ss >> px >> py >> pz >> t >> mx >> my >> mz >> polx >> poly >> polz >> wl))
0073         {
0074             G4cerr << "WARNING: skipping malformed line " << lineno << ": " << line << G4endl;
0075             continue;
0076         }
0077
0078         sphoton p = {};
0079         p.pos = {px, py, pz};
0080         p.time = t;
0081         p.mom = {mx, my, mz};
0082         p.pol = {polx, poly, polz};
0083         p.wavelength = wl;
0084         result.push_back(p);
0085     }
0086     return result;
0087 }
0088
0089 struct PrimaryGenerator : G4VUserPrimaryGeneratorAction
0090 {
0091     std::filesystem::path photon_file;
0092     SEvt *sev;
0093
0094     PrimaryGenerator(std::filesystem::path photon_file, SEvt *sev) : photon_file(photon_file), sev(sev)
0095     {
0096     }
0097
0098     void GeneratePrimaries(G4Event *event) override
0099     {
0100         std::vector<sphoton> sphotons = load_photons_txt(photon_file);
0101         if (sphotons.empty())
0102         {
0103             G4cerr << "ERROR: no photons loaded from " << photon_file << G4endl;
0104             return;
0105         }
0106
0107         G4cout << "Loaded " << sphotons.size() << " photons from " << photon_file << G4endl;
0108
0109         // Build NP array (N,4,4) for SEvt
0110         size_t num_floats = sphotons.size() * 4 * 4;
0111         float *data = reinterpret_cast<float *>(sphotons.data());
0112         NP *photons = NP::MakeFromValues<float>(data, num_floats);
0113         photons->reshape({static_cast<int64_t>(sphotons.size()), 4, 4});
0114
0115         for (const sphoton &p : sphotons)
0116         {
0117             G4ThreeVector position_mm(p.pos.x, p.pos.y, p.pos.z);
0118             G4double time_ns = p.time;
0119             G4ThreeVector direction(p.mom.x, p.mom.y, p.mom.z);
0120             G4double wavelength_nm = p.wavelength;
0121             G4ThreeVector polarization(p.pol.x, p.pol.y, p.pol.z);
0122
0123             G4PrimaryVertex *vertex = new G4PrimaryVertex(position_mm, time_ns);
0124             G4double kineticEnergy = h_Planck * c_light / (wavelength_nm * nm);
0125
0126             G4PrimaryParticle *particle = new G4PrimaryParticle(G4OpticalPhoton::Definition());
0127             particle->SetKineticEnergy(kineticEnergy);
0128             particle->SetMomentumDirection(direction);
0129             particle->SetPolarization(polarization);
0130
0131             vertex->SetPrimary(particle);
0132             event->AddPrimaryVertex(vertex);
0133         }
0134
0135         sev->SetInputPhoton(photons);
0136     }
0137 };
0138
0139 struct EventAction : G4UserEventAction
0140 {
0141     SEvt *sev;
0142     unsigned int fTotalOpticksHits{0};
0143
0144     EventAction(SEvt *sev) : sev(sev)
0145     {
0146     }
0147
0148     void BeginOfEventAction(const G4Event *event) override
0149     {
0150         sev->beginOfEvent(event->GetEventID());
0151     }
0152
0153     void EndOfEventAction(const G4Event *event) override
0154     {
0155         int eventID = event->GetEventID();
0156         sev->addEventConfigArray();
0157         sev->gather();
0158         sev->endOfEvent(eventID);
0159
0160         // GPU-based simulation
0161         G4CXOpticks *gx = G4CXOpticks::Get();
0162
0163         gx->simulate(eventID, false);
0164         cudaDeviceSynchronize();
0165
0166         SEvt *sev_gpu = SEvt::Get_EGPU();
0167         unsigned int num_hits = sev_gpu->GetNumHit(0);
0168
0169         std::cout << "Opticks: NumHits:  " << num_hits << std::endl;
0170         fTotalOpticksHits += num_hits;
0171
0172         std::ofstream outFile("opticks_hits_output.txt");
0173         if (!outFile.is_open())
0174         {
0175             std::cerr << "Error opening output file!" << std::endl;
0176         }
0177         else
0178         {
0179             for (int idx = 0; idx < int(num_hits); idx++)
0180             {
0181                 sphoton hit;
0182                 sev_gpu->getHit(hit, idx);
0183                 G4ThreeVector position = G4ThreeVector(hit.pos.x, hit.pos.y, hit.pos.z);
0184                 G4ThreeVector direction = G4ThreeVector(hit.mom.x, hit.mom.y, hit.mom.z);
0185                 G4ThreeVector polarization = G4ThreeVector(hit.pol.x, hit.pol.y, hit.pol.z);
0186                 outFile << hit.time << " " << hit.wavelength << "  " << "(" << position.x() << ", " << position.y()
0187                         << ", " << position.z() << ")  " << "(" << direction.x() << ", " << direction.y() << ", "
0188                         << direction.z() << ")  " << "(" << polarization.x() << ", " << polarization.y() << ", "
0189                         << polarization.z() << ")" << std::endl;
0190             }
0191             outFile.close();
0192         }
0193
0194         gx->reset(eventID);
0195     }
0196
0197     unsigned int GetTotalOpticksHits() const
0198     {
0199         return fTotalOpticksHits;
0200     }
0201 };
0202
0203 struct RunAction : G4UserRunAction
0204 {
0205     EventAction *fEventAction;
0206
0207     RunAction(EventAction *eventAction) : fEventAction(eventAction)
0208     {
0209     }
0210
0211     void BeginOfRunAction(const G4Run *) override
0212     {
0213     }
0214
0215     void EndOfRunAction(const G4Run *) override
0216     {
0217         std::cout << "Opticks: TotalHits:  " << fEventAction->GetTotalOpticksHits() << std::endl;
0218     }
0219 };
0220
0221 struct G4App
0222 {
0223     G4App(std::filesystem::path photon_file, std::filesystem::path gdml_file)
0224         : sev(SEvt::CreateOrReuse_ECPU()), det_cons_(new DetectorConstruction(gdml_file)),
0225           prim_gen_(new PrimaryGenerator(photon_file, sev)), event_act_(new EventAction(sev)),
0226           run_act_(new RunAction(event_act_))
0227     {
0228     }
0229
0230     SEvt *sev;
0231
0232     G4VUserDetectorConstruction *det_cons_;
0233     G4VUserPrimaryGeneratorAction *prim_gen_;
0234     EventAction *event_act_;
0235     RunAction *run_act_;
0236 };