From 24fc6e0ca4dfc3907af7df89b0e9dd1865ac1f1e Mon Sep 17 00:00:00 2001 From: mmlynari <104267016+mmlynari@users.noreply.github.com> Date: Wed, 28 Aug 2024 22:53:42 +0200 Subject: [PATCH] update CellPositionsTool for HCal Allegro (#91) * update CellPositionsTool for HCal Allegro * revert changes in CellPositionsHCalBarrelPhiThetaSegTool and add CellPositionsHCalPhiThetaSegTool * remove a typo * include comments and improve documentation * do not use GaudiAlg * replace GaudiAlg by GaudiKernel --- .../CellPositionsHCalPhiThetaSegTool.cpp | 224 ++++++++++++++++++ .../CellPositionsHCalPhiThetaSegTool.h | 83 +++++++ 2 files changed, 307 insertions(+) create mode 100644 RecFCCeeCalorimeter/src/components/CellPositionsHCalPhiThetaSegTool.cpp create mode 100644 RecFCCeeCalorimeter/src/components/CellPositionsHCalPhiThetaSegTool.h diff --git a/RecFCCeeCalorimeter/src/components/CellPositionsHCalPhiThetaSegTool.cpp b/RecFCCeeCalorimeter/src/components/CellPositionsHCalPhiThetaSegTool.cpp new file mode 100644 index 00000000..754db0f5 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CellPositionsHCalPhiThetaSegTool.cpp @@ -0,0 +1,224 @@ +#include "CellPositionsHCalPhiThetaSegTool.h" + +#include "edm4hep/CalorimeterHitCollection.h" +#include + +using dd4hep::DetElement; + +DECLARE_COMPONENT(CellPositionsHCalPhiThetaSegTool) + +CellPositionsHCalPhiThetaSegTool::CellPositionsHCalPhiThetaSegTool( + const std::string &type, + const std::string &name, + const IInterface *parent) + : AlgTool(type, name, parent) +{ + declareInterface(this); +} + +StatusCode CellPositionsHCalPhiThetaSegTool::initialize() +{ + StatusCode sc = AlgTool::initialize(); + if (sc.isFailure()) + return sc; + m_geoSvc = service("GeoSvc"); + if (!m_geoSvc) + { + error() << "Unable to locate Geometry service." << endmsg; + return StatusCode::FAILURE; + } + // get PhiTheta segmentation + m_segmentation = dynamic_cast( + 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(); + + // check if decoder contains "layer" + std::vector 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 StatusCode::FAILURE; + } + + // retrieve radii from the LayeredCalorimeterData 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(); + if (!theExtension) + { + error() << "The detector element does not have the required LayeredCalorimeterData extension." << endmsg; + return StatusCode::FAILURE; + } + + // Debug information to check the number of layers retrieved from the LayeredCalorimeterData extension + m_layersRetrieved = &(theExtension->layers); + debug() << "Number of layers retrieved: " << m_layersRetrieved->size() << endmsg; + + if (m_detectorName=="HCalBarrel"){ + m_radii = CellPositionsHCalPhiThetaSegTool::calculateLayerRadiiBarrel(); + } + else if (m_detectorName=="HCalThreePartsEndcap"){ + m_radii = CellPositionsHCalPhiThetaSegTool::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") + { + // Check that the vector containing number of layers in each cylinder is provided + if (m_numLayersHCalThreeParts.empty()) + { + error() << "The vector m_numLayersHCalThreeParts is empty." << endmsg; + return StatusCode::FAILURE; + } + // Check that the total number of layers provided in the steering file + // matches the total number of layers in the geometry xml file + for (unsigned int i=0; isize() != 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; +} + +void CellPositionsHCalPhiThetaSegTool::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 = CellPositionsHCalPhiThetaSegTool::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 << endmsg; + } + debug() << "Output positions collection size: " << outputColl.size() << endmsg; +} + + +dd4hep::Position CellPositionsHCalPhiThetaSegTool::xyzPosition(const uint64_t &aCellId) const +{ + dd4hep::DDSegmentation::CellID volumeId = aCellId; + m_decoder->set(volumeId, "phi", 0); + m_decoder->set(volumeId, "theta", 0); + + int layer = m_decoder->get(volumeId, "layer"); + // get radius in cm + double radius = m_radii[layer]; + // get local position (for r=1) + auto inSeg = m_segmentation->position(aCellId); + // scale by radius to get global position + dd4hep::Position outSeg(inSeg.x() * radius, inSeg.y() * radius, inSeg.z() * radius); + + // MM: TBD the z-coordinate still needs to be carefully validated + // at the first glance it seems to be off in some cases in the 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 CellPositionsHCalPhiThetaSegTool::layerId(const uint64_t &aCellId) +{ + int layer; + layer = m_decoder->get(aCellId, "layer"); + return layer; +} + +// calculate layer radii from LayeredCalorimeterData extension +// which is included in the geometry description +std::vector CellPositionsHCalPhiThetaSegTool::calculateLayerRadii(unsigned int startIndex, unsigned int endIndex) { + std::vector radii; + + for (unsigned int idxLayer = startIndex; idxLayer < endIndex; ++idxLayer) { + const dd4hep::rec::LayeredCalorimeterStruct::Layer & theLayer = m_layersRetrieved->at(idxLayer); + // inner radius of a given layer + double layerInnerRadius = theLayer.distance; + // radial dimension of a given layer + double layerThickness = theLayer.sensitive_thickness; + + // Calculate mid-radius for the current layer (layerInnerRadius+layerOuterRadius)/2 + double layerMidRadius = layerInnerRadius + layerThickness / 2; + + radii.push_back(layerMidRadius); + } + return radii; +} + +// calculateLayerRadii should be used for HCalTileBarrel which is formed +// by a single cylinder with a constant number of radial layers +std::vector CellPositionsHCalPhiThetaSegTool::calculateLayerRadiiBarrel() { + return calculateLayerRadii(0, m_layersRetrieved->size()); +} + +// calculateLayerRadiiEndcap should be used for HCalThreePartsEndcap which is formed +// by three cylinders along the z-coordinate and each cylinder has a different number of radial layers +std::vector CellPositionsHCalPhiThetaSegTool::calculateLayerRadiiEndcap() { + std::vector radii; + + // Calculate radii for each part (cylinder) 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 CellPositionsHCalPhiThetaSegTool::finalize() { return AlgTool::finalize(); } diff --git a/RecFCCeeCalorimeter/src/components/CellPositionsHCalPhiThetaSegTool.h b/RecFCCeeCalorimeter/src/components/CellPositionsHCalPhiThetaSegTool.h new file mode 100644 index 00000000..bc51d110 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CellPositionsHCalPhiThetaSegTool.h @@ -0,0 +1,83 @@ +#ifndef RECCALORIMETER_CellPositionsHCalPhiThetaSegTool_H +#define RECCALORIMETER_CellPositionsHCalPhiThetaSegTool_H + +// Gaudi +#include "GaudiKernel/AlgTool.h" +#include "GaudiKernel/ServiceHandle.h" + +// k4geo +#include "detectorCommon/DetUtils_k4geo.h" +#include "detectorSegmentations/FCCSWGridPhiTheta_k4geo.h" + +// k4FWCore +#include "k4FWCore/DataHandle.h" +#include "k4Interface/IGeoSvc.h" +#include "k4Interface/ICellPositionsTool.h" + +// DD4hep +#include "DD4hep/Readout.h" +#include "DD4hep/Volumes.h" +#include "DDSegmentation/Segmentation.h" +#include + +// ROOT +#include "TGeoManager.h" + +class IGeoSvc; +namespace DD4hep { +namespace DDSegmentation { +class Segmentation; +} +} + +/** @class CellPositionsHCalPhiThetaSegTool + * + * Tool to determine each Calorimeter cell position. + * + * For the FCCee HCal Barrel and EndCap with phi-theta segmentation, + * determined from the segmentation and the LayeredCalorimeterData extension. + * The LayeredCalorimeterData extension is part of the geometry description. + * + * @author Michaela Mlynarikova + */ + +class CellPositionsHCalPhiThetaSegTool : public AlgTool, virtual public ICellPositionsTool { +public: + CellPositionsHCalPhiThetaSegTool(const std::string& type, const std::string& name, const IInterface* parent); + ~CellPositionsHCalPhiThetaSegTool() = 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; + + virtual std::vector calculateLayerRadii(unsigned int startIndex, unsigned int endIndex); + + virtual std::vector calculateLayerRadiiBarrel(); + + virtual std::vector calculateLayerRadiiEndcap(); + +private: + /// Pointer to the geometry service + SmartIF m_geoSvc; + /// Name of the hadronic calorimeter readout + Gaudi::Property m_readoutName{this, "readoutName", "HCalBarrelReadout"}; + /// Name of the hadronic calorimeter + Gaudi::Property m_detectorName{this, "detectorName", "HCalBarrel"}; + /// Theta-phi segmentation + dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo* m_segmentation = nullptr; + /// Cellid decoder + dd4hep::DDSegmentation::BitFieldCoder* m_decoder = nullptr; + /// vector to store calculated layer radii + std::vector m_radii; + /// layers retrieved from the geometry + const std::vector* m_layersRetrieved = nullptr; + /// for the HCal Endcap, one needs to provide the number of layers in each cylinder + Gaudi::Property> m_numLayersHCalThreeParts{this, "numLayersHCalThreeParts", {6,9,22}}; +}; +#endif /* RECCALORIMETER_CellPositionsHCalPhiThetaSegTool_H */ \ No newline at end of file