Add tools for the implmentation of noise for ALLEGRO v3 (#107)

* Add tools for the implmentation of noise for ALLEGRO v3 * Steering file for addNoise * Conflict: remove ALLEGRO tests * Temporarily add the test for ALLEGRO addNoise (to be moved to ddsim+k4run) * Try to get the new test working --------- Co-authored-by: Zhibo Wu <[email protected]> Co-authored-by: Zhibo Wu <[email protected]>
HEP-FCC · Sep 13, 2024 · af83315 · af83315
1 parent 54c0674
commit af83315
Show file tree

Hide file tree

Showing 5 changed files with 762 additions and 4 deletions.
diff --git a/RecCalorimeter/src/components/CreateCaloCells.h b/RecCalorimeter/src/components/CreateCaloCells.h
@@ -28,7 +28,7 @@ class IGeoSvc;
 /** @class CreateCaloCells
  *
  *  Algorithm for creating calorimeter cells from Geant4 hits.
- *  Tube geometry with PhiEta segmentation expected.
+ *  Tube geometry with ModuleThetaMerged segmentation expected.
  *
  *  Flow of the program:
  *  1/ Merge Geant4 energy deposits with same cellID
@@ -70,7 +70,7 @@ class CreateCaloCells : public Gaudi::Algorithm {
   /// Handle for the calorimeter cells noise tool
   mutable ToolHandle<INoiseCaloCellsTool> m_noiseTool{"NoiseCaloCellsFlatTool", this};
   /// Handle for the geometry tool
-  ToolHandle<ICalorimeterTool> m_geoTool{"TubeLayerPhiEtaCaloTool", this};
+  ToolHandle<ICalorimeterTool> m_geoTool{"TubeLayerModuleThetaMergedCaloTool", this};
 
   /// Add crosstalk to cells?
   Gaudi::Property<bool> m_addCrosstalk{this, "addCrosstalk", false, "Add crosstalk effect?"};
@@ -92,7 +92,7 @@ class CreateCaloCells : public Gaudi::Algorithm {
   mutable DataHandle<edm4hep::CalorimeterHitCollection> m_cells{"cells", Gaudi::DataHandle::Writer, this};
   MetaDataHandle<std::string> m_cellsCellIDEncoding{m_cells, edm4hep::labels::CellIDEncoding, Gaudi::DataHandle::Writer};
   /// Name of the detector readout
-  Gaudi::Property<std::string> m_readoutName{this, "readoutName", "ECalBarrelPhiEta", "Name of the detector readout"};
+  Gaudi::Property<std::string> m_readoutName{this, "readoutName", "ECalBarrelModuleThetaMerged", "Name of the detector readout"};
   /// Name of active volumes
   Gaudi::Property<std::string> m_activeVolumeName{this, "activeVolumeName", "_sensitive", "Name of the active volumes"};
   /// Name of active layers for sampling calorimeter
@@ -116,7 +116,7 @@ class CreateCaloCells : public Gaudi::Algorithm {
    * This property won't be needed anymore.
    */
   unsigned int m_activeVolumesNumber;
-  /// Use only volume ID? If false, using PhiEtaSegmentation
+  /// Use only volume ID? If false, using ModuleThetaMergedSegmentation
   bool m_useVolumeIdOnly;
 
   /// Pointer to the geometry service

diff --git a/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.cpp b/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.cpp
@@ -0,0 +1,99 @@
+#include "TubeLayerModuleThetaMergedCaloTool.h"
+
+// segm
+#include "DD4hep/Detector.h"
+#include "DD4hep/MultiSegmentation.h"
+#include "detectorCommon/DetUtils_k4geo.h"
+#include "k4Interface/IGeoSvc.h"
+
+DECLARE_COMPONENT(TubeLayerModuleThetaMergedCaloTool)
+
+TubeLayerModuleThetaMergedCaloTool::TubeLayerModuleThetaMergedCaloTool(const std::string& type, const std::string& name,
+                                                 const IInterface* parent)
+    : AlgTool(type, name, parent), m_geoSvc("GeoSvc", name) {
+  declareInterface<ICalorimeterTool>(this);
+}
+
+StatusCode TubeLayerModuleThetaMergedCaloTool::initialize() {
+  StatusCode sc = AlgTool::initialize();
+  if (sc.isFailure()) return sc;
+
+  if (!m_geoSvc) {
+    error() << "Unable to locate Geometry Service. "
+            << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
+    return StatusCode::FAILURE;
+  }
+  if (m_readoutName != "") {
+    // Check if readouts exist
+    info() << "Readout: " << m_readoutName << endmsg;
+    if (m_geoSvc->getDetector()->readouts().find(m_readoutName) == m_geoSvc->getDetector()->readouts().end()) {
+      error() << "Readout <<" << m_readoutName << ">> does not exist." << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+  return sc;
+}
+
+StatusCode TubeLayerModuleThetaMergedCaloTool::finalize() { return AlgTool::finalize(); }
+
+StatusCode TubeLayerModuleThetaMergedCaloTool::prepareEmptyCells(std::unordered_map<uint64_t, double>& aCells) {
+  // Get the total number of active volumes in the geometry
+  auto highestVol = gGeoManager->GetTopVolume();
+  unsigned int numLayers;
+  if (!m_activeVolumesNumber) {
+    numLayers = det::utils::countPlacedVolumes(highestVol, m_activeVolumeName);
+  } else {
+    // used when MergeLayers tool is used. To be removed once MergeLayer gets replaced by RedoSegmentation.
+    numLayers = m_activeVolumesNumber;
+  }
+  info() << "Number of active layers " << numLayers << endmsg;
+
+  // get segmentation
+  dd4hep::DDSegmentation::Segmentation *aSegmentation = m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation();
+  std::string segmentationType = aSegmentation->type();
+  info() << "Segmentation type : " << segmentationType << endmsg;
+
+  dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo *moduleThetaSegmentation = nullptr;
+  if (segmentationType == "FCCSWGridModuleThetaMerged_k4geo") {
+    moduleThetaSegmentation = dynamic_cast<dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo *>(aSegmentation);
+    info() << "Segmentation: bins in Module " << moduleThetaSegmentation->nModules() << endmsg;
+  }
+  else {
+    error() << "Unable to cast segmentation pointer!!!! Tool only applicable to FCCSWGridModuleThetaMerged_k4geo segmentation." << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  // get decoder and extrema
+  auto decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder();
+  std::vector<std::pair<int, int>> extrema;
+  extrema.push_back(std::make_pair(0, numLayers - 1));
+  extrema.push_back(std::make_pair(0, 0));
+  extrema.push_back(std::make_pair(0, 0));
+
+  for (unsigned int ilayer = 0; ilayer < numLayers; ilayer++) {
+    dd4hep::DDSegmentation::CellID volumeId = 0;
+    (*decoder)["system"].set(volumeId, 4);
+    (*decoder)["layer"].set(volumeId, ilayer);
+    (*decoder)["theta"].set(volumeId, 0);
+    (*decoder)["module"].set(volumeId, 0);
+    // Get number of segmentation cells within the active volume, for given layer
+    auto numCells = det::utils::numberOfCells(volumeId, *moduleThetaSegmentation);
+    // extrema[1]: Range of module ID (0, max module ID)
+    extrema[1] = std::make_pair(0, (numCells[0] - 1) * moduleThetaSegmentation->mergedModules(ilayer));
+    // extrema[2]: Range of theta ID (min theta ID, max theta ID
+    extrema[2] = std::make_pair(numCells[2], numCells[2] + (numCells[1] - 1) * moduleThetaSegmentation->mergedThetaCells(ilayer));
+    debug() << "Layer: " << ilayer << endmsg;
+    debug() << "Number of segmentation cells in (module, theta): " << numCells << endmsg;
+    // Loop over segmentation cells
+    for (unsigned int imodule = 0; imodule < numCells[0]; imodule++) {
+      for (unsigned int itheta = 0; itheta < numCells[1]; itheta++) {
+        dd4hep::DDSegmentation::CellID cellId = volumeId; 
+	decoder->set(cellId, "module", imodule * moduleThetaSegmentation->mergedModules(ilayer));
+        decoder->set(cellId, "theta", numCells[2] + itheta * moduleThetaSegmentation->mergedThetaCells(ilayer)); // start from the minimum existing theta cell in this layer
+        aCells.insert(std::pair<dd4hep::DDSegmentation::CellID, double>(cellId, 0));
+      }
+    }
+  }
+
+  return StatusCode::SUCCESS;
+}
diff --git a/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.h b/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.h
@@ -0,0 +1,59 @@
+#ifndef RECCALORIMETER_TUBELAYERMODULETHETAMERGEDCALOTOOL_H
+#define RECCALORIMETER_TUBELAYERMODULETHETAMERGEDCALOTOOL_H
+
+// from Gaudi
+#include "GaudiKernel/AlgTool.h"
+
+// k4FWCore
+#include "k4Interface/ICalorimeterTool.h"
+
+class IGeoSvc;
+
+/** @class TubeLayerModuleThetaMergedCaloTool Reconstruction/RecCalorimeter/src/components/TubeLayerModuleThetaMergedCaloTool.h
+ *TubeLayerModuleThetaMergedCaloTool.h
+ *
+ *  Tool for geometry-dependent settings of the digitisation.
+ *  It assumes cylindrical geometry (layers) and phi-theta segmentation.
+ *
+ *  @author Anna Zaborowska
+ *  @author Zhibo Wu
+ */
+
+class TubeLayerModuleThetaMergedCaloTool : public AlgTool, virtual public ICalorimeterTool {
+public:
+  TubeLayerModuleThetaMergedCaloTool(const std::string& type, const std::string& name, const IInterface* parent);
+  virtual ~TubeLayerModuleThetaMergedCaloTool() = default;
+  virtual StatusCode initialize() final;
+  virtual StatusCode finalize() final;
+  /** Prepare a map of all existing cells in current geometry.
+   *   Active layers (cylindrical tubes) are looked in the geometry manager by name ('\b activeVolumeName').
+   *   Corresponding bitfield name is given in '\b activeFieldName'.
+   *   If users wants to limit the number of active layers, it is possible by setting '\b activeVolumesNumber'.
+   *   The total number of cells N = n_layer * n_theta * n_phi, where
+   *   n_layer is number of layers (taken from geometry if activeVolumesNumber not set),
+   *   n_theta is number of eta bins in that layer,
+   *   n_phi is number of phi bins (the same for each layer).
+   *   For more explanation please [see reconstruction documentation](@ref md_reconstruction_doc_reccalorimeter).
+   *   @param[out] aCells map of existing cells (and deposited energy, set to 0)
+   *   return Status code.
+   */
+  virtual StatusCode prepareEmptyCells(std::unordered_map<uint64_t, double>& aCells) final;
+
+private:
+  /// Pointer to the geometry service
+  ServiceHandle<IGeoSvc> m_geoSvc;
+  /// Name of the detector readout
+  Gaudi::Property<std::string> m_readoutName{this, "readoutName", ""};
+  /// Name of active volumes
+  Gaudi::Property<std::string> m_activeVolumeName{this, "activeVolumeName", "LAr_sensitive"};
+  /// Name of active layers for sampling calorimeter
+  Gaudi::Property<std::string> m_activeFieldName{this, "activeFieldName", "active_layer"};
+  /// Name of the fields describing the segmented volume
+  Gaudi::Property<std::vector<std::string>> m_fieldNames{this, "fieldNames"};
+  /// Values of the fields describing the segmented volume
+  Gaudi::Property<std::vector<int>> m_fieldValues{this, "fieldValues"};
+  /// Temporary: for use with MergeLayer tool
+  Gaudi::Property<unsigned int> m_activeVolumesNumber{this, "activeVolumesNumber", 0};
+};
+
+#endif /* RECCALORIMETER_TUBELAYERMODULETHETAMERGEDCALOTOOL_H */
diff --git a/RecFCCeeCalorimeter/CMakeLists.txt b/RecFCCeeCalorimeter/CMakeLists.txt
@@ -50,3 +50,10 @@ add_test(NAME FCCeeLAr_benchmarkCorrection
          WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
 )
 set_test_env(FCCeeLAr_benchmarkCorrection)
+
+# test for ALLEGRO addNoise, to be moved ddsim+k4run
+add_test(NAME FCCeeLAr_ecalNoise
+         COMMAND k4run RecFCCeeCalorimeter/tests/options/runEcalBarrel_addNoise_thetamodulemerged.py
+         WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+)
+set_test_env(FCCeeLAr_ecalNoise)