From 70c3ffbfdac31a2a366fade7da8c3f55b7659daa Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori Date: Thu, 21 Sep 2023 23:28:33 +0200 Subject: [PATCH] move cell position code for module-theta merged segmentation to new tool in RecFCCeeCalorimeter --- ...lPositionsECalBarrelModuleThetaSegTool.cpp | 146 ++++++++++++++++ ...ellPositionsECalBarrelModuleThetaSegTool.h | 65 +++++++ .../CellPositionsECalBarrelTool.cpp | 161 +++--------------- .../components/CellPositionsECalBarrelTool.h | 8 +- 4 files changed, 234 insertions(+), 146 deletions(-) create mode 100644 RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelModuleThetaSegTool.cpp create mode 100644 RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelModuleThetaSegTool.h diff --git a/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelModuleThetaSegTool.cpp b/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelModuleThetaSegTool.cpp new file mode 100644 index 00000000..84284c8c --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelModuleThetaSegTool.cpp @@ -0,0 +1,146 @@ +#include "CellPositionsECalBarrelModuleThetaSegTool.h" + +// EDM +#include "edm4hep/CalorimeterHitCollection.h" + +#include + +DECLARE_COMPONENT(CellPositionsECalBarrelModuleThetaSegTool) + +CellPositionsECalBarrelModuleThetaSegTool::CellPositionsECalBarrelModuleThetaSegTool(const std::string& type, const std::string& name, + const IInterface* parent) + : GaudiTool(type, name, parent) { + declareInterface(this); +} + +StatusCode CellPositionsECalBarrelModuleThetaSegTool::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 segmentation + m_segmentation = dynamic_cast(m_geoSvc->lcdd()->readout(m_readoutName).segmentation().segmentation()); + if (m_segmentation == nullptr) { + error() << "There is no module-theta segmentation!!!!" << endmsg; + return StatusCode::FAILURE; + } + debug() << "Found merged module-theta segmentation" << endmsg; + for (int iLayer=0; iLayernLayers(); iLayer++) { + info() << "Layer : " << iLayer + << " theta merge : " << m_segmentation->mergedThetaCells(iLayer) + << " module merge : " << m_segmentation->mergedModules(iLayer) << endmsg; + if (m_segmentation->mergedThetaCells(iLayer)<1) { + error() << "Number of cells merged along theta should be >= 1!!!!" << endmsg; + } + if (m_segmentation->mergedModules(iLayer)<1) { + error() << "Number of modules merged should be >= 1!!!!" << endmsg; + } + } + + + // Take readout bitfield decoder from GeoSvc + m_decoder = m_geoSvc->lcdd()->readout(m_readoutName).idSpec().decoder(); + m_volman = m_geoSvc->lcdd()->volumeManager(); + + // 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 sc; +} + +void CellPositionsECalBarrelModuleThetaSegTool::getPositions(const edm4hep::CalorimeterHitCollection& aCells, + edm4hep::CalorimeterHitCollection& outputColl) { + + debug() << "Input collection size : " << aCells.size() << endmsg; + + // Loop through input cell collection, call xyzPosition method for each cell + // and assign position to cloned hit to be saved in outputColl + for (const auto& cell : aCells) { + auto outSeg = CellPositionsECalBarrelModuleThetaSegTool::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 CellPositionsECalBarrelModuleThetaSegTool::xyzPosition(const uint64_t& aCellId) const { + + dd4hep::Position outSeg; + double radius; + + // for module-theta merged segmentation, the local position returned + // by the segmentation class is theta of group of merged cells + // and relative phi wrt first module in group of merged modules + // in a vector3 scaled such that R_xy=1 + + // find position of volume corresponding to first of group of merged cells + debug() << "cellID: " << aCellId << endmsg; + dd4hep::DDSegmentation::CellID volumeId = aCellId; + m_decoder->set(volumeId, "theta", 0); + debug() << "volumeID: " << volumeId << endmsg; + 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; + + // get R, phi of volume + radius = std::sqrt(std::pow(outGlobal[0], 2) + std::pow(outGlobal[1], 2)); + double phi = std::atan2(outGlobal[1], outGlobal[0]); + debug() << "R (mm), phi of volume : \t" << radius / dd4hep::mm << " , " << phi << endmsg; + + // now get offset in theta and in phi from cell ID (due to theta grid + merging in theta/modules) + // the local position is normalised to r_xy=1 so theta is atan(1/z) + dd4hep::DDSegmentation::Vector3D inSeg = m_segmentation->position(aCellId); + debug() << "Local position of cell (mm) : \t" + << inSeg.x() / dd4hep::mm << "\t" + << inSeg.y() / dd4hep::mm << "\t" + << inSeg.z() / dd4hep::mm << endmsg; + double dphi = std::atan2(inSeg.y(), inSeg.x()); + debug() << "Local phi of cell : \t" << dphi << endmsg; + phi += dphi; + outSeg = dd4hep::Position(radius * std::cos(phi), + radius * std::sin(phi), + radius * inSeg.z()); + debug() << "Position of cell (mm) : \t" << outSeg.x() / dd4hep::mm << "\t" << outSeg.y() / dd4hep::mm << "\t" + << outSeg.z() / dd4hep::mm << "\n" + << endmsg; + + return outSeg; +} + +int CellPositionsECalBarrelModuleThetaSegTool::layerId(const uint64_t& aCellId) { + int layer; + dd4hep::DDSegmentation::CellID cID = aCellId; + layer = m_decoder->get(cID, "layer"); + return layer; +} + +StatusCode CellPositionsECalBarrelModuleThetaSegTool::finalize() { return GaudiTool::finalize(); } diff --git a/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelModuleThetaSegTool.h b/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelModuleThetaSegTool.h new file mode 100644 index 00000000..4164a512 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelModuleThetaSegTool.h @@ -0,0 +1,65 @@ +#ifndef RECCALORIMETER_CELLPOSITIONSECALBARRELMODULETHETASEGTOOL_H +#define RECCALORIMETER_CELLPOSITIONSECALBARRELMODULETHETASEGTOOL_H + +// GAUDI +#include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/ServiceHandle.h" + +// FCCSW +#include "DetCommon/DetUtils.h" +#include "k4Interface/IGeoSvc.h" +#include "DetSegmentation/FCCSWGridModuleThetaMerged.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 CellPositionsECalBarrelModuleThetaSegTool Reconstruction/RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelModuleThetaSegTool.h + * CellPositionsECalBarrelModuleThetaSegTool.h + * + * Tool to determine each Calorimeter cell position. + * + * For the FCCee Barrel ECAL, determined from the placed volumes and the FCCSW theta-module segmentation. + * + * @author Giovanni Marchiori + */ + +class CellPositionsECalBarrelModuleThetaSegTool : public GaudiTool, virtual public ICellPositionsTool { +public: + CellPositionsECalBarrelModuleThetaSegTool(const std::string& type, const std::string& name, const IInterface* parent); + ~CellPositionsECalBarrelModuleThetaSegTool() = 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 m_geoSvc; + /// Name of the electromagnetic calorimeter readout + Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelModuleThetaMerged"}; + /// Merged module-theta segmentation + dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged* m_segmentation; + /// Cellid decoder + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + /// Volume manager + dd4hep::VolumeManager m_volman; +}; +#endif /* RECCALORIMETER_CELLPOSITIONSECALBARRELMODULETHETASEGTOOL_H */ diff --git a/RecFCChhCalorimeter/src/components/CellPositionsECalBarrelTool.cpp b/RecFCChhCalorimeter/src/components/CellPositionsECalBarrelTool.cpp index d9004d3a..67d34b49 100644 --- a/RecFCChhCalorimeter/src/components/CellPositionsECalBarrelTool.cpp +++ b/RecFCChhCalorimeter/src/components/CellPositionsECalBarrelTool.cpp @@ -3,8 +3,6 @@ // EDM #include "edm4hep/CalorimeterHitCollection.h" -#include - DECLARE_COMPONENT(CellPositionsECalBarrelTool) CellPositionsECalBarrelTool::CellPositionsECalBarrelTool(const std::string& type, const std::string& name, @@ -21,54 +19,16 @@ StatusCode CellPositionsECalBarrelTool::initialize() { error() << "Unable to locate Geometry service." << endmsg; return StatusCode::FAILURE; } - - // get segmentation - auto segmentation = m_geoSvc->lcdd()->readout(m_readoutName).segmentation().segmentation(); - m_segmentationType = -1; - - // determine segmentation type - // could probably do it using method segmentation::type()... - m_segmentation = dynamic_cast(segmentation); - if (m_segmentation != nullptr) { - info() << "Found phi-eta segmentation" << endmsg; - m_segmentationType = 0; - } - else { - warning() << "There is no phi-eta segmentation, trying phi-theta" << endmsg; - m_segmentation = dynamic_cast(segmentation); - if (m_segmentation != nullptr) { - info() << "Found phi-theta segmentation" << endmsg; - m_segmentationType = 1; - } - else { - warning() << "There is no phi-theta segmentation, trying merged module-theta merged" << endmsg; - m_segmentation = dynamic_cast(segmentation); - if (m_segmentation != nullptr) { - m_segmentationType = 2; - info() << "Found merged module-theta segmentation" << endmsg; - dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged* seg = dynamic_cast(m_segmentation); - for (int iLayer=0; iLayer<12; iLayer++) { - info() << "Layer : " << iLayer << " theta merge : " << seg->mergedThetaCells(iLayer) << " module merge : " << seg->mergedModules(iLayer) << endmsg; - if (seg->mergedThetaCells(iLayer)<1) { - error() << "Number of cells merged along theta should be >= 1!!!!" << endmsg; - } - if (seg->mergedModules(iLayer)<1) { - error() << "Number of modules merged should be >= 1!!!!" << endmsg; - } - } - } - else { - error() << "There is no module-theta segmentation!!!!" << endmsg; - return StatusCode::FAILURE; - } - } + // get PhiEta segmentation + m_segmentation = dynamic_cast( + m_geoSvc->lcdd()->readout(m_readoutName).segmentation().segmentation()); + if (m_segmentation == nullptr) { + error() << "There is no phi-eta segmentation!!!!" << endmsg; + return StatusCode::FAILURE; } - - // Take readout bitfield decoder from GeoSvc m_decoder = m_geoSvc->lcdd()->readout(m_readoutName).idSpec().decoder(); m_volman = m_geoSvc->lcdd()->volumeManager(); - // check if decoder contains "layer" std::vector fields; for (uint itField = 0; itField < m_decoder->size(); itField++) { @@ -78,7 +38,6 @@ StatusCode CellPositionsECalBarrelTool::initialize() { if (iter == fields.end()) { error() << "Readout does not contain field: 'layer'" << endmsg; } - return sc; } @@ -86,9 +45,7 @@ void CellPositionsECalBarrelTool::getPositions(const edm4hep::CalorimeterHitColl edm4hep::CalorimeterHitCollection& outputColl) { debug() << "Input collection size : " << aCells.size() << endmsg; - - // Loop through input cell collection, call xyzPosition method for each cell - // and assign position to cloned hit to be saved in outputColl + // Loop through cell collection for (const auto& cell : aCells) { auto outSeg = CellPositionsECalBarrelTool::xyzPosition(cell.getCellID()); auto edmPos = edm4hep::Vector3f(); @@ -110,97 +67,21 @@ void CellPositionsECalBarrelTool::getPositions(const edm4hep::CalorimeterHitColl } dd4hep::Position CellPositionsECalBarrelTool::xyzPosition(const uint64_t& aCellId) const { - - dd4hep::Position outSeg; double radius; - - if (m_segmentationType==2) { - - // for module-theta merged segmentation, the local position returned - // by the segmentation class is theta of group of merged cells - // and relative phi wrt first module in group of merged modules - // in a vector3 scaled such that R_xy=1 - - // find position of volume corresponding to first of group of merged cells - debug() << "cellID: " << aCellId << endmsg; - dd4hep::DDSegmentation::CellID volumeId = aCellId; - m_decoder->set(volumeId, "theta", 0); - debug() << "volumeID: " << volumeId << endmsg; - 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; - - // get R, phi of volume - radius = std::sqrt(std::pow(outGlobal[0], 2) + std::pow(outGlobal[1], 2)); - double phi = std::atan2(outGlobal[1], outGlobal[0]); - debug() << "R (mm), phi of volume : \t" << radius / dd4hep::mm << " , " << phi << endmsg; - - // now get offset in theta and in phi from cell ID (due to theta grid + merging in theta/modules) - // the local position is normalised to r_xy=1 so theta is atan(1/z) - dd4hep::DDSegmentation::Vector3D inSeg = m_segmentation->position(aCellId); - debug() << "Local position of cell (mm) : \t" - << inSeg.x() / dd4hep::mm << "\t" - << inSeg.y() / dd4hep::mm << "\t" - << inSeg.z() / dd4hep::mm << endmsg; - double dphi = std::atan2(inSeg.y(), inSeg.x()); - debug() << "Local phi of cell : \t" << dphi << endmsg; - phi += dphi; - outSeg = dd4hep::Position(radius * std::cos(phi), - radius * std::sin(phi), - radius * inSeg.z()); - /* - double dtheta = (inSeg.z() == 0.) ? M_PI_2 : std::atan(1./inSeg.z()); - if (inSeg.z() == 0.) - outSeg = dd4hep::Position(radius * std::cos(phi), - radius * std::sin(phi), - 0.0); - else - outSeg = dd4hep::Position(radius * std::cos(phi), - radius * std::sin(phi), - radius * std::cos(dtheta)/std::sin(dtheta)); - */ - } - else { - // for grid eta-phi and theta-phi segmentations, the local position returned - // by the segmentation class is actually the global position but - // scaled such that R_xy = 1 - // so first need to find the proper value of R and then to rescale the - // vector - - // find radius of volume at phi bin=0, eta/theta bin=0 - dd4hep::DDSegmentation::CellID volumeId = aCellId; - m_decoder->set(volumeId, "phi", 0); - if (m_segmentationType==0) - m_decoder->set(volumeId, "eta", 0); - else if (m_segmentationType==1) - 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 = std::sqrt(std::pow(outGlobal[0], 2) + std::pow(outGlobal[1], 2)); - - // get position of cell from segmentation class (for a radius of 1) - auto inSeg = m_segmentation->position(aCellId); - - // rescale position by actual radius - outSeg = dd4hep::Position(inSeg.x() * radius, inSeg.y() * radius, inSeg.z() * radius); - } - - debug() << "Position of cell (mm) : \t" << outSeg.x() / dd4hep::mm << "\t" << outSeg.y() / dd4hep::mm << "\t" - << outSeg.z() / dd4hep::mm << "\n" - << endmsg; + dd4hep::DDSegmentation::CellID volumeId = aCellId; + m_decoder->set(volumeId, "phi", 0); + m_decoder->set(volumeId, "eta", 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; } diff --git a/RecFCChhCalorimeter/src/components/CellPositionsECalBarrelTool.h b/RecFCChhCalorimeter/src/components/CellPositionsECalBarrelTool.h index 11f968cb..d339832e 100644 --- a/RecFCChhCalorimeter/src/components/CellPositionsECalBarrelTool.h +++ b/RecFCChhCalorimeter/src/components/CellPositionsECalBarrelTool.h @@ -9,8 +9,6 @@ #include "DetCommon/DetUtils.h" #include "k4Interface/IGeoSvc.h" #include "DetSegmentation/FCCSWGridPhiEta.h" -#include "DetSegmentation/FCCSWGridPhiTheta.h" -#include "DetSegmentation/FCCSWGridModuleThetaMerged.h" #include "k4FWCore/DataHandle.h" #include "k4Interface/ICellPositionsTool.h" @@ -57,10 +55,8 @@ class CellPositionsECalBarrelTool : public GaudiTool, virtual public ICellPositi SmartIF m_geoSvc; /// Name of the electromagnetic calorimeter readout Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelPhiEta"}; - /// Eta-phi / theta-phi / merged module-phi segmentation - dd4hep::DDSegmentation::Segmentation* m_segmentation; - /// segmentation type: -1 unknown, 0 eta-phi, 1 theta-phi, 2 theta-module merged - int m_segmentationType; + /// Eta-phi segmentation + dd4hep::DDSegmentation::FCCSWGridPhiEta* m_segmentation; /// Cellid decoder dd4hep::DDSegmentation::BitFieldCoder* m_decoder; /// Volume manager