add Calorimeter detector and output scheme

include/RMGCalorimeterDetector.hh

@@ -0,0 +1,58 @@
+// Copyright (C) 2022 Luigi Pertoldi <https://orcid.org/0000-0002-0467-2571>
+//
+// This program is free software: you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any
+// later version.
+//
+// This program is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+#ifndef _RMG_CALORIMETER_DETECTOR_HH_
+#define _RMG_CALORIMETER_DETECTOR_HH_
+
+#include <memory>
+#include <string>
+
+#include "G4VSensitiveDetector.hh"
+
+#include "RMGDetectorHit.hh"
+
+class G4Step;
+class G4HCofThisEvent;
+class G4TouchableHistory;
+/** @brief Class to describe the germanium detector, mainly handles processing of the detected hits.
+ * Extends @c G4VSensitiveDetector */
+class RMGCalorimeterDetector : public G4VSensitiveDetector {
+
+  public:
+
+    RMGCalorimeterDetector();
+    ~RMGCalorimeterDetector() = default;
+
+    RMGCalorimeterDetector(RMGCalorimeterDetector const&) = delete;
+    RMGCalorimeterDetector& operator=(RMGCalorimeterDetector const&) = delete;
+    RMGCalorimeterDetector(RMGCalorimeterDetector&&) = delete;
+    RMGCalorimeterDetector& operator=(RMGCalorimeterDetector&&) = delete;
+
+    void Initialize(G4HCofThisEvent* hit_coll) override;
+
+    /** @brief Process the detected hits computing the various parameters of a @c
+     * RMGDetectorHit and then adding this to the hit collection.*/
+    bool ProcessHits(G4Step* step, G4TouchableHistory* history) override;
+    void EndOfEvent(G4HCofThisEvent* hit_coll) override;
+
+  private:
+
+    RMGDetectorHitsCollection* fHitsCollection = nullptr;
+};
+
+
+#endif
+
+// vim: tabstop=2 shiftwidth=2 expandtab

include/RMGCalorimeterOutputScheme.hh

@@ -0,0 +1,65 @@
+// Copyright (C) 2022 Luigi Pertoldi <https://orcid.org/0000-0002-0467-2571>
+//
+// This program is free software: you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation, either version 3 of the License, or (at your option) any
+// later version.
+//
+// This program is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+#ifndef _RMG_CALORIMETER_OUTPUT_SCHEME_HH_
+#define _RMG_CALORIMETER_OUTPUT_SCHEME_HH_
+
+#include <memory>
+#include <optional>
+#include <set>
+
+#include "G4AnalysisManager.hh"
+#include "G4GenericMessenger.hh"
+
+#include "RMGCalorimeterDetector.hh"
+#include "RMGDetectorHit.hh"
+#include "RMGOutputTools.hh"
+#include "RMGVOutputScheme.hh"
+
+class G4Event;
+/** @brief Output scheme for Calorimeters.
+ *
+ *  @details This output scheme records the hits in the Calorimeters.
+ *  The properties of each @c RMGDetectorHit are recorded:
+ *  - event index,
+ *  - time,
+ *  - energy deposition,
+ */
+class RMGCalorimeterOutputScheme : public RMGVOutputScheme {
+
+  public:
+
+    RMGCalorimeterOutputScheme() {};
+
+    /** @brief Sets the names of the output columns, invoked in @c RMGRunAction::SetupAnalysisManager */
+    void AssignOutputNames(G4AnalysisManager* ana_man) override;
+
+    /** @brief Store the information from the event, invoked in @c RMGEventAction::EndOfEventAction
+     * @details Only steps with non-zero energy are stored, unless @c fDiscardZeroEnergyHits is false.
+     */
+    void StoreEvent(const G4Event* event) override;
+
+  protected:
+
+    [[nodiscard]] std::string GetNtupleNameFlat() const override { return "calorimeter"; }
+
+  private:
+
+    RMGDetectorHitsCollection* GetHitColl(const G4Event*);
+};
+
+#endif
+
+// vim: tabstop=2 shiftwidth=2 expandtab

include/RMGDetectorMetadata.hh

@@ -22,6 +22,7 @@ enum RMGDetectorType {
   kGermanium,
   kOptical,
   kScintillator,
+  kCalorimeter,
 };
 
 struct RMGDetectorMetadata {

python/remage/ipc.py

@@ -55,6 +55,7 @@
 import os
 import signal
 import subprocess
+from collections import defaultdict
 
 from ._version import __version__
 
@@ -188,7 +189,7 @@ def __init__(self, ipc_info):
         """Storage structure for the IPC messages returned by ``remage-cpp``."""
         self.ipc_info = ipc_info
 
-    def get(self, name: str, expected_len: int = 1) -> list[str]:
+    def get(self, name: str, expected_len: int = 1) -> list[str | list[str | tuple]]:
         """Return all messages of a given key ``name`` from the IPC message list.
 
         Parameters
@@ -203,9 +204,20 @@ def get(self, name: str, expected_len: int = 1) -> list[str]:
             if len(msg) == expected_len + 1 and msg[0] == name
         ]
         if expected_len == 1:
+            # remove the unneeded extra dimension
             return [msg[0] for msg in msgs]
         return msgs
 
+    def get_as_dict(self, *args, **kwargs) -> dict[str, list[str | tuple]]:
+        """Same as :meth:`.get` but return a dictionary keyed by record key."""
+        msgs = self.get(*args, **kwargs)
+
+        d = defaultdict(list)
+        for k, v in msgs:
+            d[k].append(v)
+
+        return dict(d)
+
     def get_single(self, name: str, default: str) -> str:
         """Return the single single value for the key ``name`` or ``default`` if not present.
 

python/remage/post_proc.py

@@ -1,8 +1,9 @@
 from __future__ import annotations
 
+import copy
 import logging
 import time
-from collections.abc import Sequence
+from collections.abc import Mapping, Sequence
 from contextlib import contextmanager
 from pathlib import Path
 
@@ -41,8 +42,9 @@ def post_proc(
         Path(p).suffix.lower() for p in [*remage_files, main_output_file]
     }
 
-    detector_info: list[str] = ipc_info.get("output_ntuple", 2)
-    detector_info_aux: list[str] = ipc_info.get("output_ntuple_aux", 2)
+    # these are mappings: <output scheme name> -> [<table name 1>, <table name 2>, ...]
+    detector_info = ipc_info.get_as_dict("output_ntuple", 2)
+    detector_info_aux = ipc_info.get_as_dict("output_ntuple_aux", 2)
 
     assert len(output_file_exts) == 1
 
@@ -84,19 +86,17 @@ def post_proc(
         )
         log.info(msg)
 
-        # registered scintillator or germanium detectors
-        registered_detectors = list({det[1] for det in detector_info})
-
         # extract the additional tables in the output file (not detectors)
-        extra_tables = list({det[1] for det in detector_info_aux})
 
         with tmp_renamed_files(remage_files) as original_files:
             # also get the additional tables to forward
             config = get_reboost_config(
-                registered_detectors,
-                extra_tables,
+                detector_info,
+                detector_info_aux,
                 time_window=time_window_in_us,
             )
+            msg = f"Reboost config: {config}"
+            log.debug(msg)
 
             # use reboost to post-process outputs
             reboost.build_hit(
@@ -222,20 +222,20 @@ def deduplicate_table(
 
 
 def get_reboost_config(
-    reshape_table_list: Sequence[str],
-    other_table_list: Sequence[str],
+    detector_info: Mapping[str, Sequence[str]],
+    detector_info_aux: Mapping[str, Sequence[str]],
     *,
     time_window: float = 10,
 ) -> dict:
     """Get the config file to run reboost.
 
     Parameters
     ----------
-    reshape_table_list
-        a list of the table in the remage file that need to be reshaped
-        (i.e. Germanium or Scintillator output)
-    other_table_list
-        other tables in the file.
+    detector_info
+        a mapping of tables in the remage file that need to be reshaped, keyed
+        by output scheme name (i.e. RMGGermaniumOutputScheme).
+    detector_info_aux
+        same as `detector_info`, but holds non-standard tables.
     time_window
         time window to use for building hits (in us).
 
@@ -245,24 +245,42 @@ def get_reboost_config(
     """
     config = {}
 
-    # get the config for tables to be reshaped
-    config["processing_groups"] = [
-        {
-            "name": "all",
-            "detector_mapping": [{"output": table} for table in reshape_table_list],
-            "hit_table_layout": f"reboost.shape.group.group_by_time(STEPS, {time_window})",
-            "operations": {
-                "t0": {
-                    "expression": "ak.fill_none(ak.firsts(HITS.time, axis=-1), 0)",
-                    "units": "ns",
-                },
-                "evtid": "ak.fill_none(ak.firsts(HITS.evtid, axis=-1), 0)",
-            },
-        }
+    # build a default config for all tables (exclude calorimeter tables)
+    table_list = [
+        v
+        for k, vals in detector_info.items()
+        if k != "RMGCalorimeterOutputScheme"
+        for v in vals
     ]
 
+    default_config = {
+        "name": "default",
+        "detector_mapping": [{"output": table} for table in table_list],
+        "hit_table_layout": f"reboost.shape.group.group_by_time(STEPS, {time_window})",
+        "operations": {
+            "t0": {
+                "expression": "ak.fill_none(ak.firsts(HITS.time, axis=-1), 0)",
+                "units": "ns",
+            },
+            "evtid": "ak.fill_none(ak.firsts(HITS.evtid, axis=-1), 0)",
+        },
+    }
+    config["processing_groups"] = [default_config]
+
+    if "RMGCalorimeterOutputScheme" in detector_info:
+        # special treatment for the calorimeter output scheme
+        tables = detector_info["RMGCalorimeterOutputScheme"]
+
+        calo_config = copy.deepcopy(default_config)
+        calo_config["name"] = "RMGCalorimeterOutputScheme"
+        calo_config["detector_mapping"] = [{"output": table} for table in tables]
+        # there is always one energy per event, then remove one dimension from the array
+        calo_config["operations"]["edep"] = "ak.ravel(HITS.edep)"
+
+        config["processing_groups"].append(calo_config)
+
     # forward other tables as they are
-    config["forward"] = other_table_list
+    config["forward"] = [v for vals in detector_info_aux.values() for v in vals]
 
     return config
 

src/CMakeLists.txt

@@ -5,6 +5,8 @@ set(_root ${PROJECT_SOURCE_DIR})
 
 set(PROJECT_PUBLIC_HEADERS
     ${_root}/include/RMGAnalysisReader.hh
+    ${_root}/include/RMGCalorimeterDetector.hh
+    ${_root}/include/RMGCalorimeterOutputScheme.hh
     ${_root}/include/RMGDetectorHit.hh
     ${_root}/include/RMGDetectorMetadata.hh
     ${_root}/include/RMGExceptionHandler.hh
@@ -53,6 +55,8 @@ set(PROJECT_PUBLIC_HEADERS
 
 set(PROJECT_SOURCES
     ${_root}/src/RMGAnalysisReader.cc
+    ${_root}/src/RMGCalorimeterDetector.cc
+    ${_root}/src/RMGCalorimeterOutputScheme.cc
     ${_root}/src/RMGDetectorHit.cc
     ${_root}/src/RMGExceptionHandler.cc
     ${_root}/src/RMGHardware.cc