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