examples/MC_Truth/GPUMCTruth.h

0001 #include <filesystem>
0002 #include <fstream>
0003 #include <iostream>
0004
0005 #include "G4AutoLock.hh"
0006 #include "G4BooleanSolid.hh"
0007 #include "G4Cerenkov.hh"
0008 #include "G4Electron.hh"
0009 #include "G4Event.hh"
0010 #include "G4GDMLParser.hh"
0011 #include "G4LogicalVolumeStore.hh"
0012 #include "G4OpBoundaryProcess.hh"
0013 #include "G4OpticalPhoton.hh"
0014 #include "G4PhysicalConstants.hh"
0015 #include "G4PrimaryParticle.hh"
0016 #include "G4PrimaryVertex.hh"
0017 #include "G4RunManager.hh"
0018 #include "G4RunManagerFactory.hh"
0019 #include "G4SDManager.hh"
0020 #include "G4Scintillation.hh"
0021 #include "G4SubtractionSolid.hh"
0022 #include "G4SystemOfUnits.hh"
0023 #include "G4THitsCollection.hh"
0024 #include "G4Threading.hh"
0025 #include "G4ThreeVector.hh"
0026 #include "G4Track.hh"
0027 #include "G4TrackStatus.hh"
0028 #include "G4UserEventAction.hh"
0029 #include "G4UserRunAction.hh"
0030 #include "G4UserSteppingAction.hh"
0031 #include "G4UserTrackingAction.hh"
0032 #include "G4VHit.hh"
0033 #include "G4VPhysicalVolume.hh"
0034 #include "G4VProcess.hh"
0035 #include "G4VUserDetectorConstruction.hh"
0036 #include "G4VUserPrimaryGeneratorAction.hh"
0037 #include "G4VisAttributes.hh"
0038
0039 #include "eic-opticks/g4cx/G4CXOpticks.hh"
0040 #include "eic-opticks/sysrap/NP.hh"
0041 #include "eic-opticks/sysrap/SEvt.hh"
0042 #include "eic-opticks/sysrap/STrackInfo.h"
0043 #include "eic-opticks/sysrap/spho.h"
0044 #include "eic-opticks/sysrap/sphoton.h"
0045 #include "eic-opticks/u4/U4.hh"
0046 #include "eic-opticks/u4/U4Random.hh"
0047 #include "eic-opticks/u4/U4StepPoint.hh"
0048 #include "eic-opticks/u4/U4Touchable.h"
0049 #include "eic-opticks/u4/U4Track.h"
0050
0051 namespace
0052 {
0053 G4Mutex genstep_mutex = G4MUTEX_INITIALIZER;
0054 }
0055
0056 bool IsSubtractionSolid(G4VSolid* solid)
0057 {
0058     if (!solid)
0059         return false;
0060
0061     // Check if the solid is directly a G4SubtractionSolid
0062     if (dynamic_cast<G4SubtractionSolid*>(solid))
0063         return true;
0064
0065     // If the solid is a Boolean solid, check its constituent solids
0066     G4BooleanSolid* booleanSolid = dynamic_cast<G4BooleanSolid*>(solid);
0067     if (booleanSolid)
0068     {
0069         G4VSolid* solidA = booleanSolid->GetConstituentSolid(0);
0070         G4VSolid* solidB = booleanSolid->GetConstituentSolid(1);
0071
0072         // Recursively check the constituent solids
0073         if (IsSubtractionSolid(solidA) || IsSubtractionSolid(solidB))
0074             return true;
0075     }
0076
0077     // For other solid types, return false
0078     return false;
0079 }
0080
0081 std::string str_tolower(std::string s)
0082 {
0083     std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c) { return std::tolower(c); });
0084     return s;
0085 }
0086
0087 struct PhotonHit : public G4VHit
0088 {
0089     PhotonHit() = default;
0090
0091     PhotonHit(unsigned id, G4double energy, G4double time, G4ThreeVector position, G4ThreeVector direction,
0092               G4ThreeVector polarization) :
0093         fid(id),
0094         fenergy(energy),
0095         ftime(time),
0096         fposition(position),
0097         fdirection(direction),
0098         fpolarization(polarization)
0099     {
0100     }
0101
0102     // Copy constructor
0103     PhotonHit(const PhotonHit& right) :
0104         G4VHit(right),
0105         fid(right.fid),
0106         fenergy(right.fenergy),
0107         ftime(right.ftime),
0108         fposition(right.fposition),
0109         fdirection(right.fdirection),
0110         fpolarization(right.fpolarization)
0111     {
0112     }
0113
0114     // Assignment operator
0115     const PhotonHit& operator=(const PhotonHit& right)
0116     {
0117         if (this != &right)
0118         {
0119             G4VHit::operator=(right);
0120             fid = right.fid;
0121             fenergy = right.fenergy;
0122             ftime = right.ftime;
0123             fposition = right.fposition;
0124             fdirection = right.fdirection;
0125             fpolarization = right.fpolarization;
0126         }
0127         return *this;
0128     }
0129
0130     // Equality operator
0131     G4bool operator==(const PhotonHit& right) const
0132     {
0133         return (this == &right);
0134     }
0135
0136     // Print method
0137     void Print() override
0138     {
0139         G4cout << "Detector id: " << fid << " energy: " << fenergy << " nm" << " time: " << ftime << " ns"
0140                << " position: " << fposition << " direction: " << fdirection << " polarization: " << fpolarization
0141                << G4endl;
0142     }
0143
0144     // Member variables
0145     G4int         fid{0};
0146     G4double      fenergy{0};
0147     G4double      ftime{0};
0148     G4ThreeVector fposition{0, 0, 0};
0149     G4ThreeVector fdirection{0, 0, 0};
0150     G4ThreeVector fpolarization{0, 0, 0};
0151 };
0152
0153 using PhotonHitsCollection = G4THitsCollection<PhotonHit>;
0154
0155 struct PhotonSD : public G4VSensitiveDetector
0156 {
0157     PhotonSD(G4String name) :
0158         G4VSensitiveDetector(name),
0159         fHCID(-1)
0160     {
0161         G4String HCname = name + "_HC";
0162         collectionName.insert(HCname);
0163         G4cout << collectionName.size() << "   PhotonSD name:  " << name << " collection Name: " << HCname << G4endl;
0164     }
0165
0166     void Initialize(G4HCofThisEvent* hce) override
0167     {
0168         fPhotonHitsCollection = new PhotonHitsCollection(SensitiveDetectorName, collectionName[0]);
0169         if (fHCID < 0)
0170         {
0171             // G4cout << "PhotonSD::Initialize:  " << SensitiveDetectorName << "   " << collectionName[0] << G4endl;
0172             fHCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
0173         }
0174         hce->AddHitsCollection(fHCID, fPhotonHitsCollection);
0175     }
0176
0177     G4bool ProcessHits(G4Step* aStep, G4TouchableHistory*) override
0178     {
0179         G4Track* theTrack = aStep->GetTrack();
0180         if (theTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition())
0181             return false;
0182
0183         G4double theEnergy = theTrack->GetTotalEnergy() / CLHEP::eV;
0184
0185         // Create a new hit (CopyNr is set to 0 as DetectorID is omitted)
0186         PhotonHit* newHit = new PhotonHit(
0187             0, // CopyNr set to 0
0188             theEnergy, theTrack->GetGlobalTime(), aStep->GetPostStepPoint()->GetPosition(),
0189             aStep->GetPostStepPoint()->GetMomentumDirection(), aStep->GetPostStepPoint()->GetPolarization());
0190
0191         fPhotonHitsCollection->insert(newHit);
0192         theTrack->SetTrackStatus(fStopAndKill);
0193         return true;
0194     }
0195
0196     void EndOfEvent(G4HCofThisEvent*) override
0197     {
0198         G4int NbHits = fPhotonHitsCollection->entries();
0199         G4cout << "PhotonSD::EndOfEvent Number of PhotonHits: " << NbHits << G4endl;
0200     }
0201
0202     void AddOpticksHits()
0203     {
0204         SEvt*        sev = SEvt::Get_EGPU();
0205         unsigned int num_hits = sev->GetNumHit(0);
0206
0207         for (int idx = 0; idx < int(num_hits); idx++)
0208         {
0209             sphoton hit;
0210             sev->getHit(hit, idx);
0211             G4ThreeVector position = G4ThreeVector(hit.pos.x, hit.pos.y, hit.pos.z);
0212             G4ThreeVector direction = G4ThreeVector(hit.mom.x, hit.mom.y, hit.mom.z);
0213             G4ThreeVector polarization = G4ThreeVector(hit.pol.x, hit.pol.y, hit.pol.z);
0214             int           theCreationProcessid;
0215             if (OpticksPhoton::HasCerenkovFlag(hit.flagmask))
0216             {
0217                 theCreationProcessid = 0;
0218             }
0219             else if (OpticksPhoton::HasScintillationFlag(hit.flagmask))
0220             {
0221                 theCreationProcessid = 1;
0222             }
0223             else
0224             {
0225                 theCreationProcessid = -1;
0226             }
0227             std::cout << hit.wavelength << " " << position << " " << direction << " " << polarization << std::endl;
0228
0229             PhotonHit* newHit = new PhotonHit(0, hit.wavelength, hit.time, position, direction, polarization);
0230             fPhotonHitsCollection->insert(newHit);
0231         }
0232     }
0233
0234   private:
0235     PhotonHitsCollection* fPhotonHitsCollection{nullptr};
0236     G4int                 fHCID;
0237 };
0238
0239 struct DetectorConstruction : G4VUserDetectorConstruction
0240 {
0241     DetectorConstruction(std::filesystem::path gdml_file) :
0242         gdml_file_(gdml_file)
0243     {
0244     }
0245
0246     G4VPhysicalVolume* Construct() override
0247     {
0248         parser_.Read(gdml_file_.string(), false);
0249         G4VPhysicalVolume* world = parser_.GetWorldVolume();
0250
0251         G4CXOpticks::SetGeometry(world);
0252         G4LogicalVolumeStore* lvStore = G4LogicalVolumeStore::GetInstance();
0253
0254         static G4VisAttributes invisibleVisAttr(false);
0255
0256         // Check if the store is not empty
0257         if (lvStore && !lvStore->empty())
0258         {
0259             // Iterate over all logical volumes in the store
0260             for (auto& logicalVolume : *lvStore)
0261             {
0262                 G4VSolid* solid = logicalVolume->GetSolid();
0263
0264                 // Check if the solid uses subtraction
0265                 if (IsSubtractionSolid(solid))
0266                 {
0267                     // Assign the invisible visual attributes to the logical volume
0268                     logicalVolume->SetVisAttributes(&invisibleVisAttr);
0269
0270                     // Optionally, print out the name of the logical volume
0271                     G4cout << "Hiding logical volume: " << logicalVolume->GetName() << G4endl;
0272                 }
0273             }
0274         }
0275
0276         return world;
0277     }
0278
0279     void ConstructSDandField() override
0280     {
0281         G4cout << "ConstructSDandField is called." << G4endl;
0282         G4SDManager* SDman = G4SDManager::GetSDMpointer();
0283
0284         const G4GDMLAuxMapType* auxmap = parser_.GetAuxMap();
0285         for (auto const& [logVol, listType] : *auxmap)
0286         {
0287             for (auto const& auxtype : listType)
0288             {
0289                 if (auxtype.type == "SensDet")
0290                 {
0291                     G4cout << "Attaching sensitive detector to logical volume: " << logVol->GetName() << G4endl;
0292                     G4String  name = logVol->GetName() + "_PhotonDetector";
0293                     PhotonSD* aPhotonSD = new PhotonSD(name);
0294                     SDman->AddNewDetector(aPhotonSD);
0295                     logVol->SetSensitiveDetector(aPhotonSD);
0296                 }
0297             }
0298         }
0299     }
0300
0301   private:
0302     std::filesystem::path gdml_file_;
0303     G4GDMLParser          parser_;
0304 };
0305
0306 struct PrimaryGenerator : G4VUserPrimaryGeneratorAction
0307 {
0308     SEvt* sev;
0309
0310     PrimaryGenerator(SEvt* sev) :
0311         sev(sev)
0312     {
0313     }
0314
0315     void GeneratePrimaries(G4Event* event) override
0316     {
0317         G4ThreeVector position_mm(0.0 * m, 0.0 * m, 0.0 * m);
0318         G4double      time_ns = 0;
0319         G4ThreeVector direction(0, 0.2, 0.8);
0320         G4double      wavelength_nm = 0.1;
0321
0322         G4PrimaryVertex*   vertex = new G4PrimaryVertex(position_mm, time_ns);
0323         G4PrimaryParticle* particle = new G4PrimaryParticle(G4Electron::Definition());
0324         particle->SetKineticEnergy(5 * GeV);
0325         particle->SetMomentumDirection(direction);
0326         vertex->SetPrimary(particle);
0327         event->AddPrimaryVertex(vertex);
0328     }
0329 };
0330
0331 struct EventAction : G4UserEventAction
0332 {
0333     SEvt* sev;
0334     G4int fTotalG4Hits{0};
0335
0336     EventAction(SEvt* sev) :
0337         sev(sev)
0338     {
0339     }
0340
0341     void BeginOfEventAction(const G4Event* event) override
0342     {
0343     }
0344
0345     void EndOfEventAction(const G4Event* event) override
0346     {
0347         G4HCofThisEvent* hce = event->GetHCofThisEvent();
0348         if (hce)
0349         {
0350             for (G4int i = 0; i < hce->GetNumberOfCollections(); i++)
0351             {
0352                 G4VHitsCollection* hc = hce->GetHC(i);
0353                 if (hc)
0354                 {
0355                     fTotalG4Hits += hc->GetSize();
0356                 }
0357             }
0358         }
0359     }
0360
0361     G4int GetTotalG4Hits() const
0362     {
0363         return fTotalG4Hits;
0364     }
0365 };
0366
0367 struct RunAction : G4UserRunAction
0368 {
0369     EventAction* fEventAction;
0370
0371     RunAction(EventAction* eventAction) :
0372         fEventAction(eventAction)
0373     {
0374     }
0375
0376     void BeginOfRunAction(const G4Run* run) override
0377     {
0378     }
0379
0380     void EndOfRunAction(const G4Run* run) override
0381     {
0382         if (G4Threading::IsMasterThread())
0383         {
0384             G4CXOpticks* gx = G4CXOpticks::Get();
0385
0386             auto start = std::chrono::high_resolution_clock::now();
0387             gx->simulate(0, false);
0388             cudaDeviceSynchronize();
0389             auto end = std::chrono::high_resolution_clock::now();
0390             // Compute duration
0391             std::chrono::duration<double> elapsed = end - start;
0392             std::cout << "Simulation time: " << elapsed.count() << " seconds" << std::endl;
0393
0394             // unsigned int num_hits = SEvt::GetNumHit(EGPU);
0395             SEvt*        sev = SEvt::Get_EGPU();
0396             unsigned int num_hits = sev->GetNumHit(0);
0397
0398             std::cout << "Opticks: NumCollected:  " << sev->GetNumGenstepFromGenstep(0) << std::endl;
0399             std::cout << "Opticks: NumCollected:  " << sev->GetNumPhotonCollected(0) << std::endl;
0400             std::cout << "Opticks: NumHits:  " << num_hits << std::endl;
0401             std::cout << "Geant4: NumHits:  " << fEventAction->GetTotalG4Hits() << std::endl;
0402             std::ofstream outFile("opticks_hits_output.txt");
0403             if (!outFile.is_open())
0404             {
0405                 std::cerr << "Error opening output file!" << std::endl;
0406                 return;
0407             }
0408
0409             const bool emit_trackid = getenv("OPTICKS_MC_TRUTH") != nullptr;
0410
0411             // Isolated MC-truth benchmark: two loops without file I/O so the
0412             // paired delta measures only the lookup cost (binary search +
0413             // trackid read), not the ostream formatting that dominates the
0414             // regular hit-write loop. Gated by OPTICKS_MC_TRUTH_BENCH so
0415             // normal runs pay nothing.
0416             if (getenv("OPTICKS_MC_TRUTH_BENCH"))
0417             {
0418                 volatile uint64_t sink = 0;
0419                 const NP*         hit_np = sev->getHit();
0420                 const sphoton*    hits = reinterpret_cast<const sphoton*>(hit_np->bytes());
0421                 const int         M = int(sev->gs.size());
0422                 auto              bt0 = std::chrono::high_resolution_clock::now();
0423                 for (int idx = 0; idx < int(num_hits); idx++)
0424                 {
0425                     sink ^= hits[idx].index;
0426                 }
0427                 auto bt1 = std::chrono::high_resolution_clock::now();
0428                 int  g = 0;
0429                 for (int idx = 0; idx < int(num_hits); idx++)
0430                 {
0431                     const int64_t pidx = hits[idx].index;
0432                     while (g + 1 < M && pidx >= sev->gs[g + 1].offset)
0433                         ++g;
0434                     int tid = pidx < sev->gs[g].offset + sev->gs[g].photons ? int(sev->genstep[g].trackid()) : -1;
0435                     sink ^= hits[idx].index ^ unsigned(tid);
0436                 }
0437                 auto                          bt2 = std::chrono::high_resolution_clock::now();
0438                 std::chrono::duration<double> base_t = bt1 - bt0;
0439                 std::chrono::duration<double> mc_t = bt2 - bt1;
0440                 double                        base_ns = num_hits > 0 ? 1e9 * base_t.count() / num_hits : 0.0;
0441                 double                        mc_ns = num_hits > 0 ? 1e9 * mc_t.count() / num_hits : 0.0;
0442                 std::cout << "Bench baseline:  " << base_t.count() << " s  (" << base_ns << " ns/hit)" << std::endl;
0443                 std::cout << "Bench mctruth:   " << mc_t.count() << " s  (" << mc_ns << " ns/hit)" << std::endl;
0444                 std::cout << "Bench delta:     " << (mc_t.count() - base_t.count()) << " s  (" << (mc_ns - base_ns)
0445                           << " ns/hit)   [sink=" << sink << "]" << std::endl;
0446             }
0447
0448             auto           hit_loop_start = std::chrono::high_resolution_clock::now();
0449             const NP*      hit_np_main = sev->getHit();
0450             const sphoton* hits_main = reinterpret_cast<const sphoton*>(hit_np_main->bytes());
0451             const int      M_gs = int(sev->gs.size());
0452             int            g_cur = 0;
0453             for (int idx = 0; idx < int(num_hits); idx++)
0454             {
0455                 const sphoton& hit = hits_main[idx];
0456                 G4ThreeVector  position = G4ThreeVector(hit.pos.x, hit.pos.y, hit.pos.z);
0457                 G4ThreeVector  direction = G4ThreeVector(hit.mom.x, hit.mom.y, hit.mom.z);
0458                 G4ThreeVector  polarization = G4ThreeVector(hit.pol.x, hit.pol.y, hit.pol.z);
0459                 int            theCreationProcessid;
0460                 if (OpticksPhoton::HasCerenkovFlag(hit.flagmask))
0461                 {
0462                     theCreationProcessid = 0;
0463                 }
0464                 else if (OpticksPhoton::HasScintillationFlag(hit.flagmask))
0465                 {
0466                     theCreationProcessid = 1;
0467                 }
0468                 else
0469                 {
0470                     theCreationProcessid = -1;
0471                 }
0472                 //    std::cout << "Adding hit from Opticks:" << hit.wavelength << " " << position << " " << direction
0473                 //    << "
0474                 //    "
0475                 //              << polarization << std::endl;
0476                 outFile << hit.time << " " << hit.wavelength << "  " << "(" << position.x() << ", " << position.y()
0477                         << ", " << position.z() << ")  " << "(" << direction.x() << ", " << direction.y() << ", "
0478                         << direction.z() << ")  " << "(" << polarization.x() << ", " << polarization.y() << ", "
0479                         << polarization.z() << ")  " << "CreationProcessID=" << theCreationProcessid;
0480                 if (emit_trackid)
0481                 {
0482                     const int64_t pidx = hit.index;
0483                     while (g_cur + 1 < M_gs && pidx >= sev->gs[g_cur + 1].offset)
0484                         ++g_cur;
0485                     int trackID =
0486                         pidx < sev->gs[g_cur].offset + sev->gs[g_cur].photons ? int(sev->genstep[g_cur].trackid()) : -1;
0487                     outFile << " TrackID=" << trackID;
0488                 }
0489                 outFile << std::endl;
0490             }
0491             auto                          hit_loop_end = std::chrono::high_resolution_clock::now();
0492             std::chrono::duration<double> hit_loop_elapsed = hit_loop_end - hit_loop_start;
0493             std::cout << "Hit-write loop time: " << hit_loop_elapsed.count() << " seconds"
0494                       << " (emit_trackid=" << (emit_trackid ? "1" : "0") << ", num_hits=" << num_hits << ")"
0495                       << std::endl;
0496
0497             outFile.close();
0498         }
0499     }
0500 };
0501
0502 struct SteppingAction : G4UserSteppingAction
0503 {
0504     SEvt* sev;
0505
0506     SteppingAction(SEvt* sev) :
0507         sev(sev)
0508     {
0509     }
0510
0511     void UserSteppingAction(const G4Step* aStep)
0512     {
0513         G4Track* aTrack;
0514         G4int    fNumPhotons = 0;
0515
0516         G4StepPoint*       preStep = aStep->GetPostStepPoint();
0517         G4VPhysicalVolume* volume = preStep->GetPhysicalVolume();
0518
0519         if (aStep->GetTrack()->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition())
0520         {
0521             // Kill if step count exceeds 10000 to avoid reflection forever
0522             if (aStep->GetTrack()->GetCurrentStepNumber() > 10000)
0523             {
0524                 aStep->GetTrack()->SetTrackStatus(fStopAndKill);
0525             }
0526         }
0527
0528         if (volume && volume->GetName() == "MirrorPyramid")
0529         {
0530             aTrack = aStep->GetTrack();
0531             if (aTrack->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition())
0532             {
0533                 aTrack->SetTrackStatus(fStopAndKill);
0534             }
0535         }
0536
0537         G4SteppingManager* fpSteppingManager =
0538             G4EventManager::GetEventManager()->GetTrackingManager()->GetSteppingManager();
0539         G4StepStatus stepStatus = fpSteppingManager->GetfStepStatus();
0540
0541         if (stepStatus != fAtRestDoItProc)
0542         {
0543             G4ProcessVector* procPost = fpSteppingManager->GetfPostStepDoItVector();
0544             size_t           MAXofPostStepLoops = fpSteppingManager->GetMAXofPostStepLoops();
0545             for (size_t i3 = 0; i3 < MAXofPostStepLoops; i3++)
0546             {
0547                 if ((*procPost)[i3]->GetProcessName() == "Cerenkov")
0548                 {
0549                     aTrack = aStep->GetTrack();
0550                     const G4DynamicParticle*   aParticle = aTrack->GetDynamicParticle();
0551                     G4double                   charge = aParticle->GetDefinition()->GetPDGCharge();
0552                     const G4Material*          aMaterial = aTrack->GetMaterial();
0553                     G4MaterialPropertiesTable* MPT = aMaterial->GetMaterialPropertiesTable();
0554
0555                     G4MaterialPropertyVector* Rindex = MPT->GetProperty(kRINDEX);
0556                     if (!Rindex || Rindex->GetVectorLength() == 0)
0557                     {
0558                         G4cout << "WARNING: Material has no valid RINDEX data. Skipping Cerenkov calculation."
0559                                << G4endl;
0560                         return;
0561                     }
0562
0563                     G4Cerenkov* proc = (G4Cerenkov*)(*procPost)[i3];
0564                     fNumPhotons = proc->GetNumPhotons();
0565
0566                     G4AutoLock lock(&genstep_mutex); // <-- Mutex is locked here
0567
0568                     if (fNumPhotons > 0)
0569                     {
0570                         G4double Pmin = Rindex->Energy(0);
0571                         G4double Pmax = Rindex->GetMaxEnergy();
0572                         G4double nMax = Rindex->GetMaxValue();
0573                         G4double beta1 = aStep->GetPreStepPoint()->GetBeta();
0574                         G4double beta2 = aStep->GetPostStepPoint()->GetBeta();
0575                         G4double beta = (beta1 + beta2) * 0.5;
0576                         G4double BetaInverse = 1. / beta;
0577                         G4double maxCos = BetaInverse / nMax;
0578                         G4double maxSin2 = (1.0 - maxCos) * (1.0 + maxCos);
0579                         G4double MeanNumberOfPhotons1 =
0580                             proc->GetAverageNumberOfPhotons(charge, beta1, aMaterial, Rindex);
0581                         G4double MeanNumberOfPhotons2 =
0582                             proc->GetAverageNumberOfPhotons(charge, beta2, aMaterial, Rindex);
0583
0584                         U4::CollectGenstep_G4Cerenkov_modified(aTrack, aStep, fNumPhotons, BetaInverse, Pmin, Pmax,
0585                                                                maxCos, maxSin2, MeanNumberOfPhotons1,
0586                                                                MeanNumberOfPhotons2);
0587                     }
0588                 }
0589                 if ((*procPost)[i3]->GetProcessName() == "Scintillation")
0590                 {
0591                     G4Scintillation* proc1 = (G4Scintillation*)(*procPost)[i3];
0592                     fNumPhotons = proc1->GetNumPhotons();
0593                     if (fNumPhotons > 0)
0594                     {
0595                         aTrack = aStep->GetTrack();
0596                         const G4Material*          aMaterial = aTrack->GetMaterial();
0597                         G4MaterialPropertiesTable* MPT = aMaterial->GetMaterialPropertiesTable();
0598                         if (!MPT || !MPT->ConstPropertyExists(kSCINTILLATIONTIMECONSTANT1))
0599                         {
0600                             G4cout << "WARNING: Material has no valid SCINTILLATIONTIMECONSTANT1 data. Skipping "
0601                                       "Scintillation calculation."
0602                                    << G4endl;
0603                             return;
0604                         }
0605                         G4double SCINTILLATIONTIMECONSTANT1 = MPT->GetConstProperty(kSCINTILLATIONTIMECONSTANT1);
0606
0607                         U4::CollectGenstep_DsG4Scintillation_r4695(aTrack, aStep, fNumPhotons, 1,
0608                                                                    SCINTILLATIONTIMECONSTANT1);
0609                     }
0610                 }
0611             }
0612         }
0613     }
0614 };
0615
0616 struct TrackingAction : G4UserTrackingAction
0617 {
0618     const G4Track* transient_fSuspend_track = nullptr;
0619     SEvt*          sev;
0620
0621     TrackingAction(SEvt* sev) :
0622         sev(sev)
0623     {
0624     }
0625
0626     void PreUserTrackingAction_Optical_FabricateLabel(const G4Track* track)
0627     {
0628     }
0629
0630     void PreUserTrackingAction(const G4Track* track) override
0631     {
0632     }
0633
0634     void PostUserTrackingAction(const G4Track* track) override
0635     {
0636     }
0637 };
0638
0639 struct G4App
0640 {
0641     G4App(std::filesystem::path gdml_file) :
0642         sev(SEvt::CreateOrReuse_EGPU()),
0643         det_cons_(new DetectorConstruction(gdml_file)),
0644         prim_gen_(new PrimaryGenerator(sev)),
0645         event_act_(new EventAction(sev)),
0646         run_act_(new RunAction(event_act_)),
0647         stepping_(new SteppingAction(sev)),
0648
0649         tracking_(new TrackingAction(sev))
0650     {
0651     }
0652
0653     //~G4App(){ G4CXOpticks::Finalize();}
0654
0655     // Create "global" event
0656     SEvt* sev;
0657
0658     G4VUserDetectorConstruction*   det_cons_;
0659     G4VUserPrimaryGeneratorAction* prim_gen_;
0660     EventAction*                   event_act_;
0661     RunAction*                     run_act_;
0662     SteppingAction*                stepping_;
0663     TrackingAction*                tracking_;
0664 };