Improve physics process store

examples/10.Geometries/Validate.C

@@ -16,16 +16,16 @@ Int_t Validate(const char* filename) {
 
     const auto geometryInfo = metadata->GetGeant4GeometryInfo();
 
-    geometryInfo->Print();
+    geometryInfo.Print();
 
-    if (geometryInfo->GetAllPhysicalVolumes().size() != 374) {
-        cout << "Incorrect number of physical volumes " << geometryInfo->GetAllPhysicalVolumes().size()
+    if (geometryInfo.GetAllPhysicalVolumes().size() != 374) {
+        cout << "Incorrect number of physical volumes " << geometryInfo.GetAllPhysicalVolumes().size()
              << endl;
         // assembly do not work on older geant4 versions...
         // return 1;
     }
-    if (geometryInfo->GetAllLogicalVolumes().size() != 22) {
-        cout << "Incorrect number of logical volumes " << geometryInfo->GetAllLogicalVolumes().size() << endl;
+    if (geometryInfo.GetAllLogicalVolumes().size() != 22) {
+        cout << "Incorrect number of logical volumes " << geometryInfo.GetAllLogicalVolumes().size() << endl;
         // assembly do not work on older geant4 versions...
         // return 1;
     }

include/PrimaryGeneratorAction.h

@@ -51,8 +51,6 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
     Int_t endEnergyBin;
     Double_t fSpectrumIntegral;
 
-    Int_t nBiasingVolumes;
-
     Double_t energyFactor;
 
     Double_t lastEnergy;

include/SteppingAction.h

@@ -22,10 +22,5 @@ class SteppingAction : public G4UserSteppingAction {
 
    private:
     SimulationManager* fSimulationManager;
-
-    G4String nom_vol, nom_part, nom_proc;
-    G4double ener_dep, eKin;
-
-    G4ThreeVector previousDirection;
 };
 #endif

src/Application.cxx

@@ -119,10 +119,14 @@ void Application::Run(const CommandLineParameters& commandLineParameters) {
 
     fSimulationManager->fRestGeant4Event = new TRestGeant4Event();
     fSimulationManager->fRestGeant4SubEvent = new TRestGeant4Event();
+
+    fSimulationManager->fRestGeant4Event->InitializeReferences(fSimulationManager->fRestRun);
+    fSimulationManager->fRestGeant4SubEvent->InitializeReferences(fSimulationManager->fRestRun);
+
     fSimulationManager->fRestRun->AddEventBranch(fSimulationManager->fRestGeant4SubEvent);
 
     fSimulationManager->fRestGeant4Track = new TRestGeant4Track();
-
+    fSimulationManager->fRestGeant4Track->SetEvent(fSimulationManager->fRestGeant4SubEvent);
     // choose the Random engine
     CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine);
     long seed = fSimulationManager->fRestGeant4Metadata->GetSeed();

src/DetectorConstruction.cxx

@@ -47,15 +47,14 @@ G4VPhysicalVolume* DetectorConstruction::Construct() {
 
     parser->Read(gdmlToRead, false);
 
-    auto geometryInfo = restG4Metadata->GetGeant4GeometryInfo();
-
-    geometryInfo->PopulateFromGdml(gdmlToRead);
+    restG4Metadata->fGeant4GeometryInfo.PopulateFromGdml(gdmlToRead);
 
     G4VPhysicalVolume* worldVolume = parser->GetWorldVolume();
 
-    geometryInfo->PopulateFromGeant4World(worldVolume);
+    restG4Metadata->fGeant4GeometryInfo.PopulateFromGeant4World(worldVolume);
 
-    geometryInfo->Print();
+    const auto& geometryInfo = restG4Metadata->GetGeant4GeometryInfo();
+    geometryInfo.Print();
 
     chdir(originDirectory);
 
@@ -66,9 +65,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct() {
     G4VPhysicalVolume* physicalVolume = GetPhysicalVolume(sensitiveVolume);
     if (!physicalVolume) {
         // sensitive volume was not found, perhaps the user specified a logical volume
-        auto physicalVolumes =
-            restG4Metadata->GetGeant4GeometryInfo()->GetAllPhysicalVolumesFromLogical(sensitiveVolume);
-
+        auto physicalVolumes = geometryInfo.GetAllPhysicalVolumesFromLogical(sensitiveVolume);
         if (physicalVolumes.size() == 1) {
             restG4Metadata->SetSensitiveVolume(physicalVolumes[0]);
             sensitiveVolume = (string)restG4Metadata->GetSensitiveVolume();
@@ -208,12 +205,11 @@ G4VPhysicalVolume* DetectorConstruction::Construct() {
 G4VPhysicalVolume* DetectorConstruction::GetPhysicalVolume(const G4String& physVolName) {
     G4PhysicalVolumeStore* physVolStore = G4PhysicalVolumeStore::GetInstance();
     TRestGeant4Metadata* restG4Metadata = fSimulationManager->fRestGeant4Metadata;
-
+    const auto& geometryInfo = restG4Metadata->GetGeant4GeometryInfo();
     vector<G4VPhysicalVolume*>::const_iterator physVol;
     for (physVol = physVolStore->begin(); physVol != physVolStore->end(); physVol++) {
         auto name = (*physVol)->GetName();
-        auto alternativeName =
-            (G4String)restG4Metadata->GetGeant4GeometryInfo()->GetAlternativeNameFromGeant4PhysicalName(name);
+        auto alternativeName = (G4String)geometryInfo.GetAlternativeNameFromGeant4PhysicalName(name);
 
         if (name == physVolName || alternativeName == physVolName) {
             return *physVol;

src/EventAction.cxx

@@ -250,7 +250,7 @@ void EventAction::ReOrderTrackIds(Int_t subId) {
         for (int n = 0; n < restG4Event->GetNumberOfTracks(); n++) {
             const auto& track = restG4Event->GetTrack(n);
             if (track.GetSubEventID() == subId - 1) {
-                if (track.isRadiactiveDecay()) {
+                if (track.ContainsProcess("RadioactiveDecay") /* Decay, verified on v11.0.2 */) {
                     subRestG4Event->SetSubEventTag(track.GetParticleName());
                 }
             }

src/SteppingAction.cxx

@@ -27,19 +27,30 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) {
     TRestGeant4Event* restG4Event = fSimulationManager->fRestGeant4Event;
     TRestGeant4Metadata* restG4Metadata = fSimulationManager->fRestGeant4Metadata;
 
+    const auto& geometryInfo = restG4Metadata->GetGeant4GeometryInfo();
     // Variables that describe a step are taken.
-    nom_vol = restG4Metadata->GetGeant4GeometryInfo()->GetAlternativeNameFromGeant4PhysicalName(
+    const auto& volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(
         (TString &&) aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName());
-    nom_part = aStep->GetTrack()->GetDefinition()->GetParticleName();
 
-    ener_dep = aStep->GetTotalEnergyDeposit();
-    eKin = aStep->GetTrack()->GetKineticEnergy() / keV;
+    const auto& particle = aStep->GetTrack()->GetDefinition();
+    const auto& particleID = particle->GetPDGEncoding();
+    const auto& particleName = particle->GetParticleName();
+
+    restG4Metadata->fGeant4PhysicsInfo.InsertParticleName(particleID, particleName);
+
+    const auto process = aStep->GetPostStepPoint()->GetProcessDefinedStep();
+    const auto& processID = process->GetProcessType() * 1000 + process->GetProcessSubType();
+    const auto& processName = process->GetProcessName();
+
+    restG4Metadata->fGeant4PhysicsInfo.InsertProcessName(processID, processName);
+
+    const auto totalEnergyDeposit = aStep->GetTotalEnergyDeposit() / keV;
+    const auto trackKineticEnergy = aStep->GetTrack()->GetKineticEnergy() / keV;
 
     auto sensitiveVolumeName =
-        restG4Metadata->GetGeant4GeometryInfo()->GetAlternativeNameFromGeant4PhysicalName(
-            restG4Metadata->GetSensitiveVolume());
+        geometryInfo.GetAlternativeNameFromGeant4PhysicalName(restG4Metadata->GetSensitiveVolume());
 
-    if (restTrack->GetParticleName() == "geantino" && sensitiveVolumeName.Data() == nom_vol) {
+    if (restTrack->GetParticleName() == "geantino" && sensitiveVolumeName.Data() == volumeName) {
         restG4Metadata->SetSaveAllEvents(true);
     }
 
@@ -60,51 +71,52 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) {
         exit(0);
     }
 
-    nom_proc = aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
-
     G4Track* aTrack = aStep->GetTrack();
 
     Double_t x = aTrack->GetPosition().x() / mm;
     Double_t y = aTrack->GetPosition().y() / mm;
     Double_t z = aTrack->GetPosition().z() / mm;
 
-    if ((G4String)restG4Metadata->GetSensitiveVolume() == nom_vol) {
-        restG4Event->AddEnergyToSensitiveVolume(ener_dep / keV);
+    if (restG4Metadata->GetSensitiveVolume() == volumeName) {
+        restG4Event->AddEnergyToSensitiveVolume(totalEnergyDeposit);
     }
 
-    TVector3 hitPosition(x, y, z);
-    Int_t pcsID = restTrack->GetProcessID(nom_proc);
-    Double_t hit_global_time = aStep->GetPreStepPoint()->GetGlobalTime() / second;
-    G4ThreeVector momentum = aStep->GetPreStepPoint()->GetMomentumDirection();
-    TVector3 momentumDirection = TVector3(momentum.x(), momentum.y(), momentum.z());  //.Unit();
+    const TVector3 hitPosition(x, y, z);
+    const Double_t hitGlobalTime = aStep->GetPreStepPoint()->GetGlobalTime() / second;
+    const G4ThreeVector& momentum = aStep->GetPreStepPoint()->GetMomentumDirection();
+    const TVector3 momentumDirection = TVector3(momentum.x(), momentum.y(), momentum.z());  //.Unit();
 
     Int_t volume = -1;
     Bool_t alreadyStored = false;
     // We check if the hit must be stored and keep it on restG4Track
     for (int volID = 0; volID < restG4Metadata->GetNumberOfActiveVolumes(); volID++) {
         if (restG4Event->isVolumeStored(volID)) {
-            if (restG4Metadata->GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme)
-                G4cout << "Step volume :" << nom_vol << "::("
+            if (restG4Metadata->GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) {
+                G4cout << "Step volume :" << volumeName << "::("
                        << (G4String)restG4Metadata->GetActiveVolumeName(volID) << ")" << G4endl;
+            }
 
             // We store the hit if we have activated in the config
-            Bool_t isActiveVolume = (nom_vol == (G4String)restG4Metadata->GetActiveVolumeName(volID));
-
+            const Bool_t isActiveVolume = (volumeName == restG4Metadata->GetActiveVolumeName(volID));
             if (isActiveVolume) {
                 volume = volID;
-                if (restG4Metadata->GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Extreme)
+                if (restG4Metadata->GetVerboseLevel() >=
+                    TRestStringOutput::REST_Verbose_Level::REST_Extreme) {
                     G4cout << "Storing hit" << G4endl;
-                restTrack->AddG4Hit(hitPosition, ener_dep / keV, hit_global_time, pcsID, volID, eKin,
-                                    momentumDirection);
+                }
+                restTrack->AddG4Hit(hitPosition, totalEnergyDeposit, hitGlobalTime, processID, volID,
+                                    trackKineticEnergy, momentumDirection);
                 alreadyStored = true;
+                restG4Metadata->fGeant4GeometryInfo.InsertVolumeName(volID, volumeName);
             }
         }
     }
 
     // See issue #65.
     // If the radioactive decay occurs in a non-active volume then the id will be -1
-    Bool_t isDecay = (nom_proc == (G4String) "RadioactiveDecay");
-    if (!alreadyStored && isDecay)
-        restTrack->AddG4Hit(hitPosition, ener_dep / keV, hit_global_time, pcsID, volume, eKin,
-                            momentumDirection);
-}
+    const Bool_t isDecay = (processName == (G4String) "RadioactiveDecay");
+    if (!alreadyStored && isDecay) {
+        restTrack->AddG4Hit(hitPosition, totalEnergyDeposit, hitGlobalTime, processID, volume,
+                            trackKineticEnergy, momentumDirection);
+    }
+}

test/src/examples.cxx

@@ -32,21 +32,97 @@ TEST(restG4, Example_01_NLDBD) {
     parameters.outputFile =
         thisExamplePath / "NLDBD_simulation.root";  // TODO: fix not working with local path
 
-    Application app;
-    app.Run(parameters);
-
-    // Run validation macro
-    const TString macro(thisExamplePath / "Validate.C");
-    gROOT->ProcessLine(TString::Format(".L %s", macro.Data()));  // Load macro
-    int error = 0;
-    const int result =
-        gROOT->ProcessLine(TString::Format("Validate(\"%s\")", parameters.outputFile.Data()), &error);
-    EXPECT_EQ(error, 0);
-    EXPECT_EQ(result, 0);
+    {  // Run simulation
+        Application app;
+        app.Run(parameters);
+    }
+
+    {  // Run validation macro
+        const TString macro(thisExamplePath / "Validate.C");
+        gROOT->ProcessLine(TString::Format(".L %s", macro.Data()));  // Load macro
+        int error = 0;
+        const int result =
+            gROOT->ProcessLine(TString::Format("Validate(\"%s\")", parameters.outputFile.Data()), &error);
+        EXPECT_EQ(error, 0);
+        EXPECT_EQ(result, 0);
+    }
 
     fs::current_path(originalPath);
 }
 
+TEST(restG4, TRestGeant4GeometryInfo_TRestGeant4PhysicsInfo) {
+    // Test "TRestGeant4GeometryInfo" and "TRestGeant4PhysicsInfo" even though its from Geant4Lib, we need a
+    // simulation file, so we placed the test here
+
+    const auto originalPath = fs::current_path();
+    const auto thisExamplePath = examplesPath / "01.NLDBD";
+    const auto resultsFile = thisExamplePath / "NLDBD_simulation.root";
+
+    cout << "results file: " << resultsFile << endl;
+
+    // If previous example was ran we can use generated file, otherwise we generate it again
+    if (!TRestTools::CheckFileIsAccessible(resultsFile)) {
+        cout << "Results file not found, generating file..." << endl;
+        fs::current_path(thisExamplePath);
+
+        CommandLineParameters parameters;
+        parameters.rmlFile = "NLDBD.rml";
+        parameters.outputFile = resultsFile;
+
+        Application app;
+        app.Run(parameters);
+
+        fs::current_path(originalPath);
+    } else {
+        cout << "Results file found, skipping generation" << endl;
+    }
+
+    TRestRun run(resultsFile);
+    // Test `TRestGeant4Metadata::GetUnambiguousGlobalInstance`
+    auto geant4Metadata = (TRestGeant4Metadata*)run.GetMetadataClass("TRestGeant4Metadata");
+    EXPECT_EQ(geant4Metadata != nullptr, true);
+
+    const auto& geometryInfo = geant4Metadata->GetGeant4GeometryInfo();
+
+    cout << "Printing geometry info" << endl;
+    geometryInfo.Print();
+
+    EXPECT_EQ(geometryInfo.GetAllLogicalVolumes().size() == 4, true);
+    EXPECT_EQ(geometryInfo.GetAllPhysicalVolumes().size() == 4, true);
+
+    for (const auto& logicalVolume : geometryInfo.GetAllLogicalVolumes()) {
+        if (logicalVolume == "World") {
+            const auto& physicals = geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume);
+            EXPECT_EQ(physicals.size() == 1, true);
+            const auto& physical = physicals[0];
+            EXPECT_EQ(physical == "World_PV", true);
+        } else if (logicalVolume == "gasVolume") {
+            const auto& physicals = geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume);
+            EXPECT_EQ(physicals.size() == 1, true);
+            const auto& physical = physicals[0];
+            EXPECT_EQ(physical == "gas", true);
+        } else if (logicalVolume == "vesselVolume") {
+            const auto& physicals = geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume);
+            EXPECT_EQ(physicals.size() == 1, true);
+            const auto& physical = physicals[0];
+            EXPECT_EQ(physical == "vessel", true);
+        } else if (logicalVolume == "waterTankVolume") {
+            const auto& physicals = geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume);
+            EXPECT_EQ(physicals.size() == 1, true);
+            const auto& physical = physicals[0];
+            EXPECT_EQ(physical == "waterTank", true);
+        } else {
+            cout << "unexpected logical volume found!" << endl;
+            exit(1);
+        }
+    }
+
+    const auto& physicsInfo = geant4Metadata->GetGeant4PhysicsInfo();
+
+    cout << "Printing physics info" << endl;
+    physicsInfo.Print();
+}
+
 TEST(restG4, Example_04_Muons) {
     // cd into example
     const auto originalPath = fs::current_path();
@@ -114,6 +190,16 @@ TEST(restG4, Example_07_Decay_FullChain) {
     Application app;
     app.Run(parameters);
 
+    // print processes
+    TRestRun run(parameters.outputFile.Data());
+    auto geant4Metadata = (TRestGeant4Metadata*)run.GetMetadataClass("TRestGeant4Metadata");
+    EXPECT_EQ(geant4Metadata != nullptr, true);
+
+    const auto& physicsInfo = geant4Metadata->GetGeant4PhysicsInfo();
+
+    cout << "Printing physics info" << endl;
+    physicsInfo.Print();
+
     // Run validation macro
     const TString macro(thisExamplePath / "Validate.C");
     gROOT->ProcessLine(TString::Format(".L %s", macro.Data()));  // Load macro