store detector global coordinates, deduplicate auxiliary tables and convert to LH5 structs

docs/manual/output.md

@@ -573,6 +573,41 @@ More information on how to use the TCM is provided in {ref}`manual-analysis`,
 while more documentation about how the TCM is generated is available at
 {func}`pygama.evt.tcm.generate_tcm_cols`.
 
+## Detector origins
+
+_remage_ stores the global coordinates of each Germanium detectors in a LH5
+struct called `detector_origins` (or in a table if reshaping is off):
+
+```
+/
+├── detector_origins · struct{B00000C,B00000D,...}
+│   ├── B00000C · struct{xloc,yloc,zloc}
+│   │   ├── xloc · real
+│   │   ├── yloc · real
+│   │   └── zloc · real
+│   ├── B00000D · struct{xloc,yloc,zloc}
+│   │   ├── xloc · real
+│   │   ├── yloc · real
+│   │   └── zloc · real
+│   └── ...
+└── ...
+```
+
+:::{note}
+
+For most volume types, the origin is the center of the volume, with some notable
+exceptions:
+
+- generic polycones (as used for the detectors in _legend-pygeom-hpges_) have
+  their own origin defined which is not at the center.
+- for boolean solids, the origin is the same as for the first constituent. With
+  subtractions, the origin might even be outside the volume.
+- The
+  [list in the official documentation](https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/Detector/Geometry/geomSolids.html)
+  contains all rules.
+
+:::
+
 ## The vertex table
 
 _remage_ stores data about the simulated event vertex in a table named `vtx`. In

include/RMGDetectorMetadata.hh

@@ -27,7 +27,9 @@ enum RMGDetectorType {
 struct RMGDetectorMetadata {
     RMGDetectorType type;
     int uid;
+    /** @brief name of the referenced physical volume */
     std::string name;
+    int copy_nr;
 };
 
 #endif

include/RMGGermaniumOutputScheme.hh

@@ -111,6 +111,8 @@ class RMGGermaniumOutputScheme : public RMGVOutputScheme {
       fPreClusterPars.track_energy_threshold = threshold;
     }
 
+    void EndOfRunAction(const G4Run*) override;
+
   protected:
 
     [[nodiscard]] std::string GetNtupleNameFlat() const override { return "germanium"; }
@@ -141,6 +143,8 @@ class RMGGermaniumOutputScheme : public RMGVOutputScheme {
 
     // mode of position to store
     RMGOutputTools::PositionMode fPositionMode = RMGOutputTools::PositionMode::kAverage;
+
+    std::map<std::string, G4ThreeVector> fDetectorOrigins;
 };
 
 #endif

include/RMGNavigationTools.hh

@@ -18,8 +18,11 @@
 
 #include <set>
 #include <string>
+#include <vector>
 
 #include "G4LogicalVolume.hh"
+#include "G4RotationMatrix.hh"
+#include "G4ThreeVector.hh"
 #include "G4VPhysicalVolume.hh"
 
 // TODO: write function that locates points in global coordinates by using an
@@ -85,6 +88,26 @@ namespace RMGNavigationTools {
    * number, and logs their names, copy numbers, and the names of their logical parent volumes.
    */
   void PrintListOfPhysicalVolumes();
+
+  struct VolumeTreeEntry {
+      VolumeTreeEntry() = delete;
+      VolumeTreeEntry(const VolumeTreeEntry&) = default;
+      VolumeTreeEntry(G4VPhysicalVolume* pv) { physvol = pv; }
+
+      G4VPhysicalVolume* physvol;
+
+      G4ThreeVector vol_global_translation; // origin
+      G4RotationMatrix vol_global_rotation; // identity
+      std::vector<G4RotationMatrix> partial_rotations;
+      std::vector<G4ThreeVector> partial_translations;
+  };
+
+  /**
+   * @brief find all ways to reach the world volume from a given physical volume.
+   * @param pv the physical volume to start with.
+   */
+  std::vector<VolumeTreeEntry> FindGlobalPositions(G4VPhysicalVolume* pv);
+
 } // namespace RMGNavigationTools
 
 #endif

include/RMGVertexConfinement.hh

@@ -288,19 +288,6 @@ class RMGVertexConfinement : public RMGVVertexGenerator {
 
   private:
 
-    struct VolumeTreeEntry {
-        VolumeTreeEntry() = delete;
-        VolumeTreeEntry(const VolumeTreeEntry&) = default;
-        VolumeTreeEntry(G4VPhysicalVolume* pv) { physvol = pv; }
-
-        G4VPhysicalVolume* physvol;
-
-        G4ThreeVector vol_global_translation; // origin
-        G4RotationMatrix vol_global_rotation; // identity
-        std::vector<G4RotationMatrix> partial_rotations;
-        std::vector<G4ThreeVector> partial_translations;
-    };
-
     void InitializePhysicalVolumes();
     void InitializeGeometricalVolumes(bool use_excluded_volumes);
     bool ActualGenerateVertex(G4ThreeVector& v);

python/remage/post_proc.py

@@ -7,9 +7,10 @@
 from pathlib import Path
 
 import h5py
+import numpy as np
 import pygama.evt
 import reboost
-from lgdo import lh5
+from lgdo import Array, Scalar, Struct, lh5
 from lgdo.lh5.concat import lh5concat
 
 from . import utils
@@ -70,12 +71,12 @@ def post_proc(
 
     time_start = time.time()
 
-    if not flat_output:
-        # if merging is on, write everything to a single file
-        output_files: list[str] | str = (
-            remage_files if not merge_output_files else main_output_file
-        )
+    # if merging is on, write everything to a single file
+    output_files: list[str] | str = (
+        remage_files if not merge_output_files else main_output_file
+    )
 
+    if not flat_output:
         msg = (
             "Reshaping "
             + ("and merging " if merge_output_files else "")
@@ -154,6 +155,12 @@ def post_proc(
 
         ipc_info.set("output", main_output_file)
 
+    # deduplicate entries of the process table.
+    ntuples_to_deduplicate = set(ipc_info.get("output_ntuple_deduplicate"))
+    for file in utils._to_list(output_files):
+        for table in ntuples_to_deduplicate:
+            deduplicate_table(file, table, "name", not flat_output)
+
     msg = f"Finished post-processing which took {int(time.time() - time_start)} s"
     log.info(msg)
 
@@ -186,6 +193,34 @@ def copy_links(
                         del links_group[link_name]
 
 
+def deduplicate_table(
+    file: str, table_name: str, unique_col: str, to_struct: bool
+) -> None:
+    table = lh5.read(table_name, file)
+    table_old = {
+        col: (table[col].view_as("np").copy(), table[col].attrs) for col in table
+    }
+    _, uniq_idx = np.unique(table_old[unique_col][0], return_index=True)
+    table.resize(len(uniq_idx))
+    for col, (nda, attrs) in table_old.items():
+        table[col] = Array(nda[uniq_idx], attrs=attrs)
+
+    if to_struct:
+        keys = list(set(table.keys()) - {unique_col})
+        d = {}
+        for idx in range(table.size):
+            obj = (
+                Struct({col: Scalar(table[col].nda[idx]) for col in keys})
+                if len(keys) > 1
+                else Scalar(table[keys[0]].nda[idx])
+            )
+            d[table[unique_col].nda[idx].decode("utf-8")] = obj
+
+        lh5.write(Struct(d), table_name, file, wo_mode="overwrite")
+    else:
+        lh5.write(table, table_name, file, wo_mode="overwrite")
+
+
 def get_reboost_config(
     reshape_table_list: Sequence[str],
     other_table_list: Sequence[str],

src/RMGGermaniumOutputScheme.cc

@@ -26,8 +26,10 @@
 
 #include "RMGGermaniumDetector.hh"
 #include "RMGHardware.hh"
+#include "RMGIpc.hh"
 #include "RMGLog.hh"
 #include "RMGManager.hh"
+#include "RMGNavigationTools.hh"
 #include "RMGOutputManager.hh"
 #include "RMGOutputTools.hh"
 #include "RMGTools.hh"
@@ -55,6 +57,18 @@ void RMGGermaniumOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) {
   const auto det_cons = RMGManager::Instance()->GetDetectorConstruction();
   const auto detectors = det_cons->GetDetectorMetadataMap();
 
+  auto detector_origins_id = rmg_man->CreateAndRegisterAuxNtuple(
+      "detector_origins",
+      "RMGGermaniumOutputScheme",
+      ana_man
+  );
+  ana_man->CreateNtupleSColumn(detector_origins_id, "name");
+  CreateNtupleFOrDColumn(ana_man, detector_origins_id, "xloc_in_m", fStoreSinglePrecisionPosition);
+  CreateNtupleFOrDColumn(ana_man, detector_origins_id, "yloc_in_m", fStoreSinglePrecisionPosition);
+  CreateNtupleFOrDColumn(ana_man, detector_origins_id, "zloc_in_m", fStoreSinglePrecisionPosition);
+  ana_man->FinishNtuple(detector_origins_id);
+  RMGIpc::SendIpcNonBlocking(RMGIpc::CreateMessage("output_ntuple_deduplicate", "detector_origins"));
+
   std::set<int> registered_uids;
   std::map<std::string, int> registered_ntuples;
   for (auto&& det : detectors) {
@@ -64,6 +78,17 @@ void RMGGermaniumOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) {
     auto had_uid = registered_uids.emplace(det.second.uid);
     if (!had_uid.second) continue;
 
+    auto pv = RMGNavigationTools::FindPhysicalVolume(det.second.name, det.second.copy_nr);
+    auto trees = RMGNavigationTools::FindGlobalPositions(pv);
+    if (trees.size() > 1) {
+      RMGLog::Out(
+          RMGLog::fatal,
+          "more than one way to reach world volume from detector ",
+          det.second.name
+      );
+    }
+    fDetectorOrigins.insert({det.second.name, trees[0].vol_global_translation});
+
     auto ntuple_name = this->GetNtupleName(det.second);
     auto ntuple_reg = registered_ntuples.find(ntuple_name);
     if (ntuple_reg != registered_ntuples.end()) {
@@ -340,6 +365,25 @@ std::optional<bool> RMGGermaniumOutputScheme::StackingActionNewStage(const int s
   return ShouldDiscardEvent(event) ? std::make_optional(false) : std::nullopt;
 }
 
+void RMGGermaniumOutputScheme::EndOfRunAction(const G4Run*) {
+  auto rmg_man = RMGOutputManager::Instance();
+  if (!rmg_man->IsPersistencyEnabled() ||
+      (G4Threading::IsMasterThread() && !RMGManager::Instance()->IsExecSequential()))
+    return;
+
+  const auto ana_man = G4AnalysisManager::Instance();
+  auto ntuple_id = rmg_man->GetAuxNtupleID("detector_origins");
+
+  for (const auto& [det, v] : fDetectorOrigins) {
+    int col_id = 0;
+    ana_man->FillNtupleSColumn(ntuple_id, col_id++, det);
+    FillNtupleFOrDColumn(ana_man, ntuple_id, col_id++, v.getX() / u::m, fStoreSinglePrecisionPosition);
+    FillNtupleFOrDColumn(ana_man, ntuple_id, col_id++, v.getY() / u::m, fStoreSinglePrecisionPosition);
+    FillNtupleFOrDColumn(ana_man, ntuple_id, col_id++, v.getZ() / u::m, fStoreSinglePrecisionPosition);
+    ana_man->AddNtupleRow(ntuple_id);
+  }
+}
+
 void RMGGermaniumOutputScheme::SetPositionModeString(std::string mode) {
 
   try {

src/RMGHardware.cc

@@ -300,7 +300,7 @@ void RMGHardware::RegisterDetector(
   fActiveDetectors.insert(type);
 
   // FIXME: can this be done with emplace?
-  auto r_value = fDetectorMetadata.insert({{pv_name, copy_nr}, {type, uid, pv_name}});
+  auto r_value = fDetectorMetadata.insert({{pv_name, copy_nr}, {type, uid, pv_name, copy_nr}});
   if (!r_value.second) { // if insertion did not take place
     RMGLog::OutFormat(
         RMGLog::warning,

src/RMGNavigationTools.cc

@@ -16,12 +16,14 @@
 #include "RMGNavigationTools.hh"
 
 #include <map>
+#include <queue>
 #include <set>
 
 #include "G4LogicalVolume.hh"
 #include "G4LogicalVolumeStore.hh"
 #include "G4Material.hh"
 #include "G4PhysicalVolumeStore.hh"
+#include "G4TransportationManager.hh"
 #include "G4UnitsTable.hh"
 
 #include "RMGLog.hh"
@@ -157,4 +159,77 @@ void RMGNavigationTools::PrintListOfPhysicalVolumes() {
   RMGLog::Out(RMGLog::summary, "Total: ", volumes.size(), " volumes");
 }
 
+
+std::vector<RMGNavigationTools::VolumeTreeEntry> RMGNavigationTools::FindGlobalPositions(
+    G4VPhysicalVolume* pv
+) {
+  auto world_volume = G4TransportationManager::GetTransportationManager()
+                          ->GetNavigatorForTracking()
+                          ->GetWorldVolume();
+
+  std::vector<VolumeTreeEntry> trees;
+  // queue for paths to the mother volume that still have to be searched.
+  std::queue<VolumeTreeEntry> q;
+  q.emplace(pv);
+
+  for (; !q.empty(); q.pop()) {
+    auto v = q.front();
+
+    if (!v.physvol)
+      RMGLog::OutDev(
+          RMGLog::fatal,
+          "nullptr detected in loop condition, this is unexpected. ",
+          "Blame RMGNavigationTools::FindDirectMother?"
+      );
+
+    v.partial_rotations.push_back(v.physvol->GetObjectRotationValue());
+    v.partial_translations.push_back(v.physvol->GetObjectTranslation());
+
+    v.vol_global_rotation = v.partial_rotations.back() * v.vol_global_rotation;
+
+    for (auto m : RMGNavigationTools::FindDirectMothers(v.physvol)) {
+      if (m != world_volume) {
+        auto v_m = VolumeTreeEntry(v); // create a copy of the current helper object.
+        v_m.physvol = m;
+        q.push(v_m);
+      } else { // we finished that branch!
+        trees.push_back(v);
+      }
+    }
+  }
+
+  RMGLog::OutFormatDev(
+      RMGLog::debug,
+      "Found {} ways to reach world volume from {}",
+      trees.size(),
+      pv->GetName()
+  );
+
+  // finalize all found paths to the mother volume.
+  for (auto&& v : trees) {
+    // world volume not included in loop
+    v.partial_translations.emplace_back(); // origin
+    v.partial_rotations.emplace_back();    // identity
+
+    // partial_rotations[0] and partial_translations[0] refer to the target
+    // volume partial_rotations[1] and partial_translations[1], to the direct
+    // mother, etc. It is necessary to rotate with respect to the frame of the
+    // mother. If there are no rotations (or only the target volume is
+    // rotated): rotations are identity matrices and vol_global_translation =
+    // sum(partial_translations)
+    for (size_t i = 0; i < v.partial_translations.size() - 1; i++) {
+      G4ThreeVector tmp = v.partial_translations[i];
+      for (size_t j = i + 1; j < v.partial_rotations.size() - 1; j++) {
+        tmp *= v.partial_rotations[j];
+      }
+      v.vol_global_translation += tmp;
+    }
+  }
+
+  if (trees.empty())
+    RMGLog::OutDev(RMGLog::fatal, "No path to world volume found, that should not be!");
+
+  return trees;
+}
+
 // vim: tabstop=2 shiftwidth=2 expandtab

src/RMGTrackOutputScheme.cc

@@ -20,6 +20,7 @@
 #include "G4EventManager.hh"
 #include "G4OpticalPhoton.hh"
 
+#include "RMGIpc.hh"
 #include "RMGLog.hh"
 #include "RMGOutputManager.hh"
 
@@ -53,6 +54,7 @@ void RMGTrackOutputScheme::AssignOutputNames(G4AnalysisManager* ana_man) {
   ana_man->CreateNtupleIColumn(pid, "procid");
   ana_man->CreateNtupleSColumn(pid, "name");
   ana_man->FinishNtuple(pid);
+  RMGIpc::SendIpcNonBlocking(RMGIpc::CreateMessage("output_ntuple_deduplicate", "processes"));
 }
 
 void RMGTrackOutputScheme::TrackingActionPre(const G4Track* track) {