Add distance to nearest surface field to GermaniumOutputScheme (#182)

* add distance to nearest surface field to GermaniumOutputScheme * Add new output field to tests
legend-exp · Dec 13, 2024 · c670c7d · c670c7d
1 parent 2202460
commit c670c7d
Show file tree

Hide file tree

Showing 9 changed files with 60 additions and 7 deletions.
diff --git a/include/RMGGermaniumDetector.hh b/include/RMGGermaniumDetector.hh
@@ -48,6 +48,7 @@ class RMGGermaniumDetectorHit : public G4VHit {
     int detector_uid = -1;
     int particle_type = -1;
     double energy_deposition = -1;
+    double distance_to_surface = -1;
     G4ThreeVector global_position;
     double global_time = -1;
 };

diff --git a/src/RMGGermaniumDetector.cc b/src/RMGGermaniumDetector.cc
@@ -19,6 +19,7 @@
 #include <stdexcept>
 #include <string>
 
+#include "G4AffineTransform.hh"
 #include "G4Circle.hh"
 #include "G4HCofThisEvent.hh"
 #include "G4OpticalPhoton.hh"
@@ -89,10 +90,14 @@ bool RMGGermaniumDetector::ProcessHits(G4Step* step, G4TouchableHistory* /*histo
 
   // we're going to use info from the pre-step point
   const auto prestep = step->GetPreStepPoint();
+  const auto position = prestep->GetPosition();
 
   // locate us
-  const auto pv_name = prestep->GetTouchableHandle()->GetVolume()->GetName();
+  const auto pv = prestep->GetTouchableHandle()->GetVolume();
+  const auto pv_name = pv->GetName();
   const auto pv_copynr = prestep->GetTouchableHandle()->GetCopyNumber();
+  const auto lv = pv->GetLogicalVolume();
+  const auto sv = lv->GetSolid();
 
   // check if physical volume is registered as germanium detector
   const auto det_cons = RMGManager::Instance()->GetDetectorConstruction();
@@ -119,9 +124,33 @@ bool RMGGermaniumDetector::ProcessHits(G4Step* step, G4TouchableHistory* /*histo
   hit->detector_uid = det_uid;
   hit->particle_type = step->GetTrack()->GetDefinition()->GetPDGEncoding();
   hit->energy_deposition = step->GetTotalEnergyDeposit();
-  hit->global_position = prestep->GetPosition();
+  hit->global_position = position;
   hit->global_time = prestep->GetGlobalTime();
 
+  // Get distance to surface.
+  // Check distance to surfaces of Mother volume
+  //
+
+  // First transform coordinates into local system
+  G4AffineTransform tf(pv->GetRotation(), pv->GetTranslation());
+  tf.Invert();
+  double dist = sv->DistanceToOut(tf.TransformPoint(position));
+
+  // Also check distance to daughters if there are any. Analogue to G4NormalNavigation.cc
+  int local_no_daughters = lv->GetNoDaughters();
+  // sample_no has to be signed, so auto typing does not work
+  for (int sample_no = local_no_daughters - 1; sample_no >= 0; sample_no--) {
+    const auto sample_physical = lv->GetDaughter(sample_no);
+    G4AffineTransform sample_tf(sample_physical->GetRotation(), sample_physical->GetTranslation());
+    sample_tf.Invert();
+    const auto sample_point = sample_tf.TransformPoint(position);
+    const auto sample_solid = sample_physical->GetLogicalVolume()->GetSolid();
+    const double sample_dist = sample_solid->DistanceToIn(sample_point);
+    if (sample_dist < dist) { dist = sample_dist; }
+  }
+
+  hit->distance_to_surface = dist;
+
   // register the hit in the hit collection for the event
   fHitsCollection->insert(hit);
 

diff --git a/src/RMGGermaniumOutputScheme.cc b/src/RMGGermaniumOutputScheme.cc
@@ -69,6 +69,7 @@ void RMGGermaniumOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) {
     CreateNtupleFOrDColumn(ana_man, id, "xloc_in_m", fStoreSinglePrecisionPosition);
     CreateNtupleFOrDColumn(ana_man, id, "yloc_in_m", fStoreSinglePrecisionPosition);
     CreateNtupleFOrDColumn(ana_man, id, "zloc_in_m", fStoreSinglePrecisionPosition);
+    CreateNtupleFOrDColumn(ana_man, id, "dist_to_surf_in_m", fStoreSinglePrecisionPosition);
 
     ana_man->FinishNtuple(id);
   }
@@ -159,6 +160,8 @@ void RMGGermaniumOutputScheme::StoreEvent(const G4Event* event) {
           fStoreSinglePrecisionPosition);
       FillNtupleFOrDColumn(ana_man, ntupleid, col_id++, hit->global_position.getZ() / u::m,
           fStoreSinglePrecisionPosition);
+      FillNtupleFOrDColumn(ana_man, ntupleid, col_id++, hit->distance_to_surface / u::m,
+          fStoreSinglePrecisionPosition);
 
       // NOTE: must be called here for hit-oriented output
       ana_man->AddNtupleRow(ntupleid);

diff --git a/tests/output/macros/ntuple-flat.hdf5ls b/tests/output/macros/ntuple-flat.hdf5ls
@@ -6,6 +6,9 @@ stp/germanium/columns
 stp/germanium/det_uid
 stp/germanium/det_uid/entries
 stp/germanium/det_uid/pages
+stp/germanium/dist_to_surf_in_m
+stp/germanium/dist_to_surf_in_m/entries
+stp/germanium/dist_to_surf_in_m/pages
 stp/germanium/edep_in_keV
 stp/germanium/edep_in_keV/entries
 stp/germanium/edep_in_keV/pages

diff --git a/tests/output/macros/ntuple-flat.hdf5ls-lh5 b/tests/output/macros/ntuple-flat.hdf5ls-lh5
@@ -1,6 +1,7 @@
 stp : struct{germanium,scintillator,vertices}
-stp/germanium : table{evtid,det_uid,particle,edep,time,xloc,yloc,zloc}
+stp/germanium : table{evtid,det_uid,particle,edep,time,xloc,yloc,zloc,dist_to_surf}
 stp/germanium/det_uid : array<1>{real}
+stp/germanium/dist_to_surf : array<1>{real}
 stp/germanium/edep : array<1>{real}
 stp/germanium/evtid : array<1>{real}
 stp/germanium/particle : array<1>{real}

diff --git a/tests/output/macros/ntuple-per-det-vol.hdf5ls b/tests/output/macros/ntuple-per-det-vol.hdf5ls
@@ -3,6 +3,9 @@ header
 stp
 stp/det1
 stp/det1/columns
+stp/det1/dist_to_surf_in_m
+stp/det1/dist_to_surf_in_m/entries
+stp/det1/dist_to_surf_in_m/pages
 stp/det1/edep_in_keV
 stp/det1/edep_in_keV/entries
 stp/det1/edep_in_keV/pages
@@ -29,6 +32,9 @@ stp/det1/zloc_in_m/entries
 stp/det1/zloc_in_m/pages
 stp/det2
 stp/det2/columns
+stp/det2/dist_to_surf_in_m
+stp/det2/dist_to_surf_in_m/entries
+stp/det2/dist_to_surf_in_m/pages
 stp/det2/edep_in_keV
 stp/det2/edep_in_keV/entries
 stp/det2/edep_in_keV/pages

diff --git a/tests/output/macros/ntuple-per-det-vol.hdf5ls-lh5 b/tests/output/macros/ntuple-per-det-vol.hdf5ls-lh5
@@ -1,13 +1,15 @@
 stp : struct{det1,det2,scint1,scint2,vertices}
-stp/det1 : table{evtid,particle,edep,time,xloc,yloc,zloc}
+stp/det1 : table{evtid,particle,edep,time,xloc,yloc,zloc,dist_to_surf}
+stp/det1/dist_to_surf : array<1>{real}
 stp/det1/edep : array<1>{real}
 stp/det1/evtid : array<1>{real}
 stp/det1/particle : array<1>{real}
 stp/det1/time : array<1>{real}
 stp/det1/xloc : array<1>{real}
 stp/det1/yloc : array<1>{real}
 stp/det1/zloc : array<1>{real}
-stp/det2 : table{evtid,particle,edep,time,xloc,yloc,zloc}
+stp/det2 : table{evtid,particle,edep,time,xloc,yloc,zloc,dist_to_surf}
+stp/det2/dist_to_surf : array<1>{real}
 stp/det2/edep : array<1>{real}
 stp/det2/evtid : array<1>{real}
 stp/det2/particle : array<1>{real}

diff --git a/tests/output/macros/ntuple-per-det.hdf5ls b/tests/output/macros/ntuple-per-det.hdf5ls
@@ -85,6 +85,9 @@ stp/det002/zloc_pre_in_m/entries
 stp/det002/zloc_pre_in_m/pages
 stp/det011
 stp/det011/columns
+stp/det011/dist_to_surf_in_m
+stp/det011/dist_to_surf_in_m/entries
+stp/det011/dist_to_surf_in_m/pages
 stp/det011/edep_in_keV
 stp/det011/edep_in_keV/entries
 stp/det011/edep_in_keV/pages
@@ -111,6 +114,9 @@ stp/det011/zloc_in_m/entries
 stp/det011/zloc_in_m/pages
 stp/det012
 stp/det012/columns
+stp/det012/dist_to_surf_in_m
+stp/det012/dist_to_surf_in_m/entries
+stp/det012/dist_to_surf_in_m/pages
 stp/det012/edep_in_keV
 stp/det012/edep_in_keV/entries
 stp/det012/edep_in_keV/pages

diff --git a/tests/output/macros/ntuple-per-det.hdf5ls-lh5 b/tests/output/macros/ntuple-per-det.hdf5ls-lh5
@@ -25,15 +25,17 @@ stp/det002/yloc_post : array<1>{real}
 stp/det002/yloc_pre : array<1>{real}
 stp/det002/zloc_post : array<1>{real}
 stp/det002/zloc_pre : array<1>{real}
-stp/det011 : table{evtid,particle,edep,time,xloc,yloc,zloc}
+stp/det011 : table{evtid,particle,edep,time,xloc,yloc,zloc,dist_to_surf}
+stp/det011/dist_to_surf : array<1>{real}
 stp/det011/edep : array<1>{real}
 stp/det011/evtid : array<1>{real}
 stp/det011/particle : array<1>{real}
 stp/det011/time : array<1>{real}
 stp/det011/xloc : array<1>{real}
 stp/det011/yloc : array<1>{real}
 stp/det011/zloc : array<1>{real}
-stp/det012 : table{evtid,particle,edep,time,xloc,yloc,zloc}
+stp/det012 : table{evtid,particle,edep,time,xloc,yloc,zloc,dist_to_surf}
+stp/det012/dist_to_surf : array<1>{real}
 stp/det012/edep : array<1>{real}
 stp/det012/evtid : array<1>{real}
 stp/det012/particle : array<1>{real}