update CellPositionsTool for HCal Allegro

RecFCCeeCalorimeter/src/components/CellPositionsHCalBarrelPhiThetaSegTool.cpp

@@ -1,6 +1,10 @@
 #include "CellPositionsHCalBarrelPhiThetaSegTool.h"
 
 #include "edm4hep/CalorimeterHitCollection.h"
+#include <DDRec/DetectorData.h>
+
+using dd4hep::DetElement;
+
 
 DECLARE_COMPONENT(CellPositionsHCalBarrelPhiThetaSegTool)
 
@@ -35,6 +39,7 @@ StatusCode CellPositionsHCalBarrelPhiThetaSegTool::initialize()
   // 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++)
@@ -46,6 +51,60 @@ StatusCode CellPositionsHCalBarrelPhiThetaSegTool::initialize()
   {
     error() << "Readout does not contain field: 'layer'" << endmsg;
   }
+
+  // needed to retrieve radii from the dd4hep extension
+  dd4hep::Detector* detector = m_geoSvc->getDetector();
+  if (!detector)
+  {
+    error() << "Unable to retrieve the detector." << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  DetElement caloDetElem = detector->detector(m_detectorName);
+  if (!caloDetElem.isValid())
+  {
+    error() << "Unable to retrieve the detector element: " << m_detectorName << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  dd4hep::rec::LayeredCalorimeterData* theExtension = caloDetElem.extension<dd4hep::rec::LayeredCalorimeterData>();
+  if (!theExtension)
+  {
+    error() << "The detector element does not have the required LayeredCalorimeterData extension." << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  // Debug information to verify the layers retrieved
+  // m_layersRetrieved is a pointer, then it needs to point to the address of the vector
+  m_layersRetrieved = &(theExtension->layers);
+  debug() << "Number of layers retrieved: " << m_layersRetrieved->size() << endmsg;
+
+  if (m_detectorName=="HCalBarrel"){
+    m_radii = CellPositionsHCalBarrelPhiThetaSegTool::calculateLayerRadiiBarrel();
+  }
+  else if (m_detectorName=="HCalThreePartsEndcap"){
+    m_radii = CellPositionsHCalBarrelPhiThetaSegTool::calculateLayerRadiiEndcap();
+  }
+  else{
+    error() << "Provided detector name in m_detectorName " << m_detectorName << " is not matching, expected inputs are HCalBarrel or HCalThreePartsEndcap" << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  unsigned int numLayersProvided = 0;
+
+  if (m_detectorName=="HCalThreePartsEndcap")
+  {
+    for (unsigned int i=0; i<m_numLayersHCalThreeParts.size(); i++) 
+    {
+      numLayersProvided += m_numLayersHCalThreeParts[i];
+    }
+
+    if (m_layersRetrieved->size() != numLayersProvided)
+    {
+      error() << "Total number of radial layers provided in m_numLayersHCalThreeParts " << numLayersProvided << " does not match the numbers retrieved from the xml file " << m_layersRetrieved->size() << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
   return sc;
 }
 
@@ -75,6 +134,7 @@ void CellPositionsHCalBarrelPhiThetaSegTool::getPositions(const edm4hep::Calorim
   debug() << "Output positions collection size: " << outputColl.size() << endmsg;
 }
 
+
 dd4hep::Position CellPositionsHCalBarrelPhiThetaSegTool::xyzPosition(const uint64_t &aCellId) const
 {
   dd4hep::DDSegmentation::CellID volumeId = aCellId;
@@ -84,22 +144,76 @@ dd4hep::Position CellPositionsHCalBarrelPhiThetaSegTool::xyzPosition(const uint6
   int layer = m_decoder->get(volumeId, "layer");
   // get radius in cm
   double radius = m_radii[layer];
-  debug() << "Layer : " << layer << "\tradius : " << radius << " cm" << endmsg;
   // get local position (for r=1)
   auto inSeg = m_segmentation->position(aCellId);
-  debug() << "Local position : x = " << inSeg.x() << " y = "  << inSeg.y() << " z = "  << inSeg.z() << endmsg;
-
   // scale by radius to get global position
   dd4hep::Position outSeg(inSeg.x() * radius, inSeg.y() * radius, inSeg.z() * radius);
+
+  // MM: the z-coordinate is not correct for Endcap 
+  debug() << "Layer : " << layer << "\tradius : " << radius << " cm" << endmsg;
+  debug() << "Local position : x = " << inSeg.x() << " y = "  << inSeg.y() << " z = "  << inSeg.z() << endmsg;
   debug() << "Global position : x = " << outSeg.x() << " y = "  << outSeg.y() << " z = "  << outSeg.z() << endmsg;
+
   return outSeg;
 }
 
+
 int CellPositionsHCalBarrelPhiThetaSegTool::layerId(const uint64_t &aCellId)
 {
   int layer;
   layer = m_decoder->get(aCellId, "layer");
   return layer;
 }
 
+std::vector<double> CellPositionsHCalBarrelPhiThetaSegTool::calculateLayerRadii(unsigned int startIndex, unsigned int endIndex) {
+  double totalDepth = 0;
+  std::vector<double> radii;
+
+  for (unsigned int idxLayer = startIndex; idxLayer < endIndex; ++idxLayer) {
+    const dd4hep::rec::LayeredCalorimeterStruct::Layer & theLayer = m_layersRetrieved->at(idxLayer);
+    double distanceFirstLayer = theLayer.distance;
+    double layerThickness = theLayer.sensitive_thickness;
+
+    // Initialize totalDepth for the first layer in the range
+    if (idxLayer == startIndex) {
+      totalDepth = distanceFirstLayer + layerThickness;
+    } else {
+      totalDepth += layerThickness;
+    }
+
+    // Calculate radius for the current layer
+    double radiusLayer = totalDepth - layerThickness / 2;
+    radii.push_back(radiusLayer);
+  }
+  return radii;
+}
+
+// to be used for HCalTileBarrel 
+std::vector<double> CellPositionsHCalBarrelPhiThetaSegTool::calculateLayerRadiiBarrel() {
+  return calculateLayerRadii(0, m_layersRetrieved->size());
+}
+
+// to be used for HCalThreePartsEndcap
+std::vector<double> CellPositionsHCalBarrelPhiThetaSegTool::calculateLayerRadiiEndcap() {
+  std::vector<double> radii;
+
+  // Calculate radii for each part and merge the results
+  unsigned int upperIndexLayerRadiiPartOne = m_numLayersHCalThreeParts[0];
+  unsigned int upperIndexLayerRadiiPartTwo = upperIndexLayerRadiiPartOne + m_numLayersHCalThreeParts[1];
+
+  // Part 1
+  auto partOneRadii = calculateLayerRadii(0, upperIndexLayerRadiiPartOne);
+  radii.insert(radii.end(), partOneRadii.begin(), partOneRadii.end());
+
+  // Part 2
+  auto partTwoRadii = calculateLayerRadii(upperIndexLayerRadiiPartOne, upperIndexLayerRadiiPartTwo);
+  radii.insert(radii.end(), partTwoRadii.begin(), partTwoRadii.end());
+
+  // Part 3
+  auto partThreeRadii = calculateLayerRadii(upperIndexLayerRadiiPartTwo, m_layersRetrieved->size());
+  radii.insert(radii.end(), partThreeRadii.begin(), partThreeRadii.end());
+
+  return radii;
+}
+
 StatusCode CellPositionsHCalBarrelPhiThetaSegTool::finalize() { return GaudiTool::finalize(); }

RecFCCeeCalorimeter/src/components/CellPositionsHCalBarrelPhiThetaSegTool.h

@@ -18,6 +18,7 @@
 #include "DD4hep/Readout.h"
 #include "DD4hep/Volumes.h"
 #include "DDSegmentation/Segmentation.h"
+#include <DDRec/DetectorData.h>
 
 // ROOT
 #include "TGeoManager.h"
@@ -48,16 +49,24 @@ class CellPositionsHCalBarrelPhiThetaSegTool : public GaudiTool, virtual public
 
   virtual int layerId(const uint64_t& aCellId) final;
 
+  virtual std::vector<double> calculateLayerRadii(unsigned int startIndex, unsigned int endIndex); 
+  virtual std::vector<double> calculateLayerRadiiBarrel(); 
+  virtual std::vector<double> calculateLayerRadiiEndcap(); 
+
 private:
   /// Pointer to the geometry service
   SmartIF<IGeoSvc> m_geoSvc;
-  /// Name of the electromagnetic calorimeter readout
+  /// Name of the calorimeter readout
   Gaudi::Property<std::string> m_readoutName{this, "readoutName", "HCalBarrelReadout"};
-  /// segmentation for FCCee cell dimensions: 4x50mm, 6x100mm, 3x200mm
-  /// GM: risky to put the numbers here... better to get them from the geometry in memory...
-  Gaudi::Property<std::vector<double>> m_radii{this, "radii", {283.55, 288.55, 293.55, 298.55, 306.05, 316.05, 326.05, 336.05, 346.05, 356.05, 371.05, 391.05, 411.05}};
+  /// Name of the calorimeter 
+  Gaudi::Property<std::string> m_detectorName{this, "detectorName", "HCalBarrel"};
+
+  std::vector<double> m_radii;
   dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo* m_segmentation;
   dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
   dd4hep::VolumeManager m_volman;
+  // layer radii calculated on the flight from the geometry 
+  const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer>* m_layersRetrieved;
+  Gaudi::Property<std::vector<int>> m_numLayersHCalThreeParts{this, "numLayersHCalThreeParts", {6,9,22}};
 };
-#endif /* RECCALORIMETER_CellPositionsHCalBarrelPhiThetaSegTool_H */
+#endif /* RECCALORIMETER_CellPositionsHCalBarrelPhiThetaSegTool_H */