Combined ECAL+HCAL topoclustering (#61)

* make warning message more explicit * add HCal barrel to tool creating noise level map * add code to create neighbours for HCAL and to link ECAL-HCAL barrels * cell positioning tool for possible ECal barrel with phi-theta segmentation
HEP-FCC · Feb 10, 2024 · 08e9f33 · 08e9f33
1 parent 01e00e6
commit 08e9f33
Show file tree

Hide file tree

Showing 7 changed files with 733 additions and 248 deletions.
diff --git a/RecCalorimeter/src/components/ConstNoiseTool.cpp b/RecCalorimeter/src/components/ConstNoiseTool.cpp
@@ -34,12 +34,12 @@ StatusCode ConstNoiseTool::initialize() {
   m_systemNoiseConstMap.emplace(11, m_HCalFwdThreshold );
 
   m_systemNoiseOffsetMap.emplace(5, 0. );
-  m_systemNoiseOffsetMap.emplace(6, 0.  );
+  m_systemNoiseOffsetMap.emplace(6, 0. );
   m_systemNoiseOffsetMap.emplace(7, 0. );
-  m_systemNoiseOffsetMap.emplace(8, 0.  );
-  m_systemNoiseOffsetMap.emplace(9, 0.  );
-  m_systemNoiseOffsetMap.emplace(10,0.  );
-  m_systemNoiseOffsetMap.emplace(11,0.  );
+  m_systemNoiseOffsetMap.emplace(8, 0. );
+  m_systemNoiseOffsetMap.emplace(9, 0. );
+  m_systemNoiseOffsetMap.emplace(10,0. );
+  m_systemNoiseOffsetMap.emplace(11,0. );
 
   // Get GeoSvc
   m_geoSvc = service("GeoSvc");
@@ -67,10 +67,10 @@ double ConstNoiseTool::getNoiseConstantPerCell(uint64_t aCellId) {
   dd4hep::DDSegmentation::CellID cID = aCellId;
   unsigned cellSystem = m_decoder->get(cID, "system");
   // cell noise in system
-  if (m_systemNoiseConstMap[cellSystem])
+  if (m_systemNoiseConstMap.count(cellSystem))
     Noise = m_systemNoiseConstMap[cellSystem];
   else
-    warning() << "No noise constants set for this subsystem! Noise of cell set to 0. " << endmsg;
+    warning() << "No noise constant set for subsystem " << cellSystem << "! Noise constant of cell set to 0. " << endmsg;
   return Noise;
 }
 
@@ -81,9 +81,9 @@ double ConstNoiseTool::getNoiseOffsetPerCell(uint64_t aCellId) {
   dd4hep::DDSegmentation::CellID cID = aCellId;
   unsigned cellSystem = m_decoder->get(cID, "system");
   // cell noise in system
-  if (m_systemNoiseOffsetMap[cellSystem])
+  if (m_systemNoiseOffsetMap.count(cellSystem))
     Noise = m_systemNoiseOffsetMap[cellSystem];
   else
-    warning() << "No noise constants set for this subsystem! Noise of cell set to 0. " << endmsg;
+    warning() << "No noise offset set for subsystem " << cellSystem << "! Noise offset of cell set to 0. " << endmsg;
   return Noise;
 }
diff --git a/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelPhiThetaSegTool.cpp b/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelPhiThetaSegTool.cpp
@@ -0,0 +1,96 @@
+#include "CellPositionsECalBarrelPhiThetaSegTool.h"
+
+// EDM
+#include "edm4hep/CalorimeterHitCollection.h"
+
+DECLARE_COMPONENT(CellPositionsECalBarrelPhiThetaSegTool)
+
+CellPositionsECalBarrelPhiThetaSegTool::CellPositionsECalBarrelPhiThetaSegTool(const std::string& type, const std::string& name,
+                                                         const IInterface* parent)
+    : GaudiTool(type, name, parent) {
+  declareInterface<ICellPositionsTool>(this);
+}
+
+StatusCode CellPositionsECalBarrelPhiThetaSegTool::initialize() {
+  StatusCode sc = GaudiTool::initialize();
+  if (sc.isFailure()) return sc;
+  m_geoSvc = service("GeoSvc");
+  if (!m_geoSvc) {
+    error() << "Unable to locate Geometry service." << endmsg;
+    return StatusCode::FAILURE;
+  }
+  // get phi-theta segmentation
+  m_segmentation = dynamic_cast<dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo*>(
+      m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation());
+  if (m_segmentation == nullptr) {
+    error() << "There is no phi-theta segmentation!!!!" << endmsg;
+    return StatusCode::FAILURE;
+  }
+  // Take readout bitfield decoder from GeoSvc
+  m_decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder();
+  m_volman = m_geoSvc->getDetector()->volumeManager();
+  // check if decoder contains "layer"
+  std::vector<std::string> fields;
+  for (uint itField = 0; itField < m_decoder->size(); itField++) {
+    fields.push_back((*m_decoder)[itField].name());
+  }
+  auto iter = std::find(fields.begin(), fields.end(), "layer");
+  if (iter == fields.end()) {
+    error() << "Readout does not contain field: 'layer'" << endmsg;
+  }
+  return sc;
+}
+
+void CellPositionsECalBarrelPhiThetaSegTool::getPositions(const edm4hep::CalorimeterHitCollection& aCells,
+                                               edm4hep::CalorimeterHitCollection& outputColl) {
+
+  debug() << "Input collection size : " << aCells.size() << endmsg;
+  // Loop through cell collection
+  for (const auto& cell : aCells) {
+    auto outSeg = CellPositionsECalBarrelPhiThetaSegTool::xyzPosition(cell.getCellID());
+    auto edmPos = edm4hep::Vector3f();
+    edmPos.x = outSeg.x() / dd4hep::mm;
+    edmPos.y = outSeg.y() / dd4hep::mm;
+    edmPos.z = outSeg.z() / dd4hep::mm;
+
+    auto positionedHit = cell.clone();
+    positionedHit.setPosition(edmPos);
+    outputColl.push_back(positionedHit);
+
+    // Debug information about cell position
+    debug() << "Cell energy (GeV) : " << positionedHit.getEnergy() << "\tcellID " << positionedHit.getCellID() << endmsg;
+    debug() << "Position of cell (mm) : \t" << outSeg.x() / dd4hep::mm << "\t" << outSeg.y() / dd4hep::mm << "\t"
+            << outSeg.z() / dd4hep::mm << "\n"
+            << endmsg;
+  }
+  debug() << "Output positions collection size: " << outputColl.size() << endmsg;
+}
+
+dd4hep::Position CellPositionsECalBarrelPhiThetaSegTool::xyzPosition(const uint64_t& aCellId) const {
+  double radius;
+  dd4hep::DDSegmentation::CellID volumeId = aCellId;
+  m_decoder->set(volumeId, "phi", 0);
+  m_decoder->set(volumeId, "theta", 0);
+  auto detelement = m_volman.lookupDetElement(volumeId);
+  const auto& transformMatrix = detelement.nominal().worldTransformation();
+  double outGlobal[3];
+  double inLocal[] = {0, 0, 0};
+  transformMatrix.LocalToMaster(inLocal, outGlobal);
+  //debug() << "Position of volume (mm) : \t" << outGlobal[0] / dd4hep::mm << "\t" << outGlobal[1] / dd4hep::mm << "\t"
+  //        << outGlobal[2] / dd4hep::mm << endmsg;
+  // radius calculated from segmenation + z postion of volumes
+  auto inSeg = m_segmentation->position(aCellId);
+  radius = std::sqrt(std::pow(outGlobal[0], 2) + std::pow(outGlobal[1], 2));
+  dd4hep::Position outSeg(inSeg.x() * radius, inSeg.y() * radius, inSeg.z() * radius);
+
+  return outSeg;
+}
+
+int CellPositionsECalBarrelPhiThetaSegTool::layerId(const uint64_t& aCellId) {
+  int layer;
+  dd4hep::DDSegmentation::CellID cID = aCellId;
+  layer = m_decoder->get(cID, "layer");
+  return layer;
+}
+
+StatusCode CellPositionsECalBarrelPhiThetaSegTool::finalize() { return GaudiTool::finalize(); }
diff --git a/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelPhiThetaSegTool.h b/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelPhiThetaSegTool.h
@@ -0,0 +1,65 @@
+#ifndef RECCALORIMETER_CELLPOSITIONSECALBARRELPHITHETASEGTOOL_H
+#define RECCALORIMETER_CELLPOSITIONSECALBARRELPHITHETASEGTOOL_H
+
+// GAUDI
+#include "GaudiAlg/GaudiTool.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+// FCCSW
+#include "detectorCommon/DetUtils_k4geo.h"
+#include "k4Interface/IGeoSvc.h"
+#include "detectorSegmentations/FCCSWGridPhiTheta_k4geo.h"
+#include "k4FWCore/DataHandle.h"
+#include "k4Interface/ICellPositionsTool.h"
+
+// DD4hep
+#include "DD4hep/Readout.h"
+#include "DD4hep/Volumes.h"
+#include "DDSegmentation/Segmentation.h"
+#include "TGeoManager.h"
+
+class IGeoSvc;
+namespace DD4hep {
+namespace DDSegmentation {
+class Segmentation;
+}
+}
+
+/** @class CellPositionsECalBarrelPhiThetaSegTool k4RecCalorimeter/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelPhiThetaSegTool.h
+ * CellPositionsECalBarrelPhiThetaSegTool.h
+ *
+ *  Tool to determine each Calorimeter cell position.
+ *
+ *   For the FCCee Barrel ECAL with phi-theta segmentation, determined from the placed volumes and the segmentation.   
+ * 
+ *  @author Giovanni Marchiori
+ */
+
+class CellPositionsECalBarrelPhiThetaSegTool : public GaudiTool, virtual public ICellPositionsTool {
+public:
+  CellPositionsECalBarrelPhiThetaSegTool(const std::string& type, const std::string& name, const IInterface* parent);
+  ~CellPositionsECalBarrelPhiThetaSegTool() = default;
+
+  virtual StatusCode initialize() final;
+
+  virtual StatusCode finalize() final;
+
+  virtual void getPositions(const edm4hep::CalorimeterHitCollection& aCells, edm4hep::CalorimeterHitCollection& outputColl) final;
+
+  virtual dd4hep::Position xyzPosition(const uint64_t& aCellId) const final;
+
+  virtual int layerId(const uint64_t& aCellId) final;
+
+private:
+  /// Pointer to the geometry service
+  SmartIF<IGeoSvc> m_geoSvc;
+  /// Name of the electromagnetic calorimeter readout
+  Gaudi::Property<std::string> m_readoutName{this, "readoutName", "ECalBarrelPhiTheta"};
+  /// Eta-phi segmentation
+  dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo* m_segmentation;
+  /// Cellid decoder
+  dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+  /// Volume manager
+  dd4hep::VolumeManager m_volman;
+};
+#endif /* RECCALORIMETER_CELLPOSITIONSECALBARRELPHITHETASEGTOOL_H */