Merge pull request #132 from rest-for-physics/print-by-volume

Method to print only some volumes
rest-for-physics · Oct 9, 2024 · e16e5cd · e16e5cd
2 parents bb6fd19 + 4f5682a
commit e16e5cd
Show file tree

Hide file tree

Showing 4 changed files with 93 additions and 4 deletions.
diff --git a/inc/TRestGeant4Event.h b/inc/TRestGeant4Event.h
@@ -242,6 +242,7 @@ class TRestGeant4Event : public TRestEvent {
     /// maxTracks : number of tracks to print, 0 = all
     void PrintActiveVolumes() const;
     void PrintEvent(int maxTracks = 0, int maxHits = 0) const;
+    void PrintEventFilterVolumes(const std::set<std::string>& volumeNames) const;
 
     inline TPad* DrawEvent(const TString& option = "") override { return DrawEvent(option, true); }
     TPad* DrawEvent(const TString& option, Bool_t autoBoundaries);

diff --git a/inc/TRestGeant4Track.h b/inc/TRestGeant4Track.h
@@ -124,6 +124,7 @@ class TRestGeant4Track {
 
     /// Prints the track information. N number of hits to print, 0 = all
     void PrintTrack(size_t maxHits = 0) const;
+    void PrintTrackFilterVolumes(const std::set<std::string>& filterVolumes) const;
 
     inline void RemoveHits() { fHits.RemoveHits(); }
 

diff --git a/src/TRestGeant4Event.cxx b/src/TRestGeant4Event.cxx
@@ -1166,6 +1166,33 @@ void TRestGeant4Event::PrintEvent(int maxTracks, int maxHits) const {
     }
 }
 
+void TRestGeant4Event::PrintEventFilterVolumes(const std::set<std::string>& volumeNames) const {
+    TRestEvent::PrintEvent();
+
+    cout << "- Total deposited energy: " << ToEnergyString(fTotalDepositedEnergy) << endl;
+    cout << "- Sensitive detectors total energy: " << ToEnergyString(fSensitiveVolumeEnergy) << endl;
+
+    cout << "- Primary source position: " << VectorToString(fPrimaryPosition) << " mm" << endl;
+
+    for (unsigned int i = 0; i < GetNumberOfPrimaries(); i++) {
+        const char* sourceNumberString =
+            GetNumberOfPrimaries() <= 1 ? ""
+                                        : TString::Format(" %d of %zu", i + 1, GetNumberOfPrimaries()).Data();
+        cout << "   - Source" << sourceNumberString << " primary particle: " << fPrimaryParticleNames[i]
+             << endl;
+        cout << "     Direction: " << VectorToString(fPrimaryDirections[i]) << endl;
+        cout << "     Energy: " << ToEnergyString(fPrimaryEnergies[i]) << endl;
+    }
+
+    cout << endl;
+    cout << "Total number of tracks: " << GetNumberOfTracks() << endl;
+
+    int nTracks = GetNumberOfTracks();
+    for (int i = 0; i < nTracks; i++) {
+        GetTrack(i).PrintTrackFilterVolumes(volumeNames);
+    }
+}
+
 Bool_t TRestGeant4Event::ContainsProcessInVolume(Int_t processID, Int_t volumeID) const {
     for (const auto& track : fTracks) {
         if (track.ContainsProcessInVolume(processID, volumeID)) {

diff --git a/src/TRestGeant4Track.cxx b/src/TRestGeant4Track.cxx
@@ -118,10 +118,11 @@ void TRestGeant4Track::PrintTrack(size_t maxHits) const {
         << (GetParentTrack() != nullptr
                 ? TString::Format(" - Parent particle: %s", GetParentTrack()->GetParticleName().Data()).Data()
                 : "")
-        << " - Created by '" << fCreatorProcess
-        << "' with initial "
-           "KE of "
-        << ToEnergyString(fInitialKineticEnergy) << "" << endl;
+        << " - Created by '" << fCreatorProcess << "' in volume '" << GetInitialVolume()
+        << "' with initial KE of " << ToEnergyString(fInitialKineticEnergy) << " - Initial position "
+        << VectorToString(fInitialPosition) << " mm at time " << ToTimeString(fGlobalTimestamp)
+        << " - Time length of " << ToTimeString(fTimeLength) << " and spatial length of "
+        << ToLengthString(fLength) << endl;
 
     cout << "   Initial position " << VectorToString(fInitialPosition) << " mm at time "
          << ToTimeString(fGlobalTimestamp) << " - Time length of " << ToTimeString(fTimeLength)
@@ -157,6 +158,65 @@ void TRestGeant4Track::PrintTrack(size_t maxHits) const {
     }
 }
 
+void TRestGeant4Track::PrintTrackFilterVolumes(const std::set<std::string>& volumeNames) const {
+    const TRestGeant4Metadata* metadata = GetGeant4Metadata();
+    if (metadata == nullptr) {
+        return;
+    }
+    bool skip = true;
+    for (unsigned int i = 0; i < GetNumberOfHits(); i++) {
+        // check volumeName is in set
+        TString volumeName = metadata->GetGeant4GeometryInfo().GetVolumeFromID(fHits.GetHitVolume(i));
+        if (volumeName.IsNull()) {
+            // in case process name is not found, use ID
+            volumeName = TString(std::to_string(fHits.GetHitVolume(i)));
+        }
+        if (volumeNames.find(volumeName.Data()) != volumeNames.end()) {
+            skip = false;
+            break;
+        }
+    }
+
+    if (skip) {
+        return;
+    }
+
+    cout
+        << " * TrackID: " << fTrackID << " - Particle: " << fParticleName << " - ParentID: " << fParentID
+        << ""
+        << (GetParentTrack() != nullptr
+                ? TString::Format(" - Parent particle: %s", GetParentTrack()->GetParticleName().Data()).Data()
+                : "")
+        << " - Created by '" << fCreatorProcess << "' in volume '" << GetInitialVolume()
+        << "' with initial KE of " << ToEnergyString(fInitialKineticEnergy) << " - Initial position "
+        << VectorToString(fInitialPosition) << " mm at time " << ToTimeString(fGlobalTimestamp)
+        << " - Time length of " << ToTimeString(fTimeLength) << " and spatial length of "
+        << ToLengthString(fLength) << endl;
+
+    size_t nHits = GetNumberOfHits();
+    for (unsigned int i = 0; i < nHits; i++) {
+        TString processName = GetProcessName(fHits.GetHitProcess(i));
+        if (processName.IsNull()) {
+            processName =
+                TString(std::to_string(fHits.GetHitProcess(i)));  // in case process name is not found, use ID
+        }
+
+        TString volumeName = metadata->GetGeant4GeometryInfo().GetVolumeFromID(fHits.GetHitVolume(i));
+        if (volumeName.IsNull()) {
+            // in case process name is not found, use ID
+            volumeName = TString(std::to_string(fHits.GetHitVolume(i)));
+        }
+        if (volumeNames.find(volumeName.Data()) == volumeNames.end()) {
+            continue;
+        }
+        cout << "      - Hit " << i << " - Energy: " << ToEnergyString(fHits.GetEnergy(i))
+             << " - Process: " << processName << " - Volume: " << volumeName
+             << " - Position: " << VectorToString(TVector3(fHits.GetX(i), fHits.GetY(i), fHits.GetZ(i)))
+             << " mm - Time: " << ToTimeString(fHits.GetTime(i))
+             << " - KE: " << ToEnergyString(fHits.GetKineticEnergy(i)) << endl;
+    }
+}
+
 Bool_t TRestGeant4Track::ContainsProcessInVolume(Int_t processID, Int_t volumeID) const {
     for (unsigned int i = 0; i < GetNumberOfHits(); i++) {
         if (fHits.GetHitProcess(i) != processID) continue;