From 08e9f336dbdf8be94c453225cdb1d881268d2811 Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori <39376142+giovannimarchiori@users.noreply.github.com> Date: Sat, 10 Feb 2024 21:15:39 +0100 Subject: [PATCH] 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 --- .../src/components/ConstNoiseTool.cpp | 18 +- ...CellPositionsECalBarrelPhiThetaSegTool.cpp | 96 ++++ .../CellPositionsECalBarrelPhiThetaSegTool.h | 65 +++ .../components/CreateFCCeeCaloNeighbours.cpp | 506 ++++++++++++------ .../components/CreateFCCeeCaloNeighbours.h | 32 +- .../CreateFCCeeCaloNoiseLevelMap.cpp | 252 ++++++--- .../components/CreateFCCeeCaloNoiseLevelMap.h | 12 +- 7 files changed, 733 insertions(+), 248 deletions(-) create mode 100644 RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelPhiThetaSegTool.cpp create mode 100644 RecFCCeeCalorimeter/src/components/CellPositionsECalBarrelPhiThetaSegTool.h diff --git a/RecCalorimeter/src/components/ConstNoiseTool.cpp b/RecCalorimeter/src/components/ConstNoiseTool.cpp index 7c64d9f1..dc007dde 100644 --- 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 new file mode 100644 index 00000000..94d222dc --- /dev/null +++ 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(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( + 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 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 new file mode 100644 index 00000000..9a4b3987 --- /dev/null +++ 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 m_geoSvc; + /// Name of the electromagnetic calorimeter readout + Gaudi::Property 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 */ diff --git a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.cpp b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.cpp index 26273200..71f9bf1a 100644 --- a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.cpp +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.cpp @@ -48,23 +48,20 @@ StatusCode CreateFCCeeCaloNeighbours::initialize() std::pair extremaECalLastLayerPhi; std::pair extremaECalLastLayerModule; std::pair extremaECalLastLayerTheta; - // double eCalThetaOffset = 0; - // double eCalThetaSize = 0; - // double eCalPhiOffset = 0; - // double eCalModuleSize = 0; - // double hCalThetaOffset = 0; - // double hCalThetaSize = 0; - // double hCalPhiOffset = 0; + int ecalBarrelId = -1; // index of ecal barrel in segmentation list dd4hep::DDSegmentation::BitFieldCoder *decoderECalBarrel = nullptr; + dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo *hcalBarrelSegmentation = nullptr; // will be used for volume connecting std::pair extremaHCalFirstLayerPhi; std::pair extremaHCalFirstLayerTheta; - std::pair extremaHCalFirstLayerZ; + int hcalBarrelId = -1; // index of hcal barrel in segmentation list dd4hep::DDSegmentation::BitFieldCoder *decoderHCalBarrel = nullptr; + dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo *ecalBarrelModuleThetaSegmentation = nullptr; + dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo *ecalBarrelPhiThetaSegmentation = nullptr; for (uint iSys = 0; iSys < m_readoutNamesSegmented.size(); iSys++) { - // Check if readouts exist + // Check if readout exists info() << "Readout: " << m_readoutNamesSegmented[iSys] << endmsg; if (m_geoSvc->getDetector()->readouts().find(m_readoutNamesSegmented[iSys]) == m_geoSvc->getDetector()->readouts().end()) { @@ -72,7 +69,7 @@ StatusCode CreateFCCeeCaloNeighbours::initialize() return StatusCode::FAILURE; } - // get Theta-based segmentation + // get theta-based segmentation dd4hep::DDSegmentation::Segmentation *aSegmentation = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).segmentation().segmentation(); if (aSegmentation == nullptr) { @@ -122,31 +119,46 @@ StatusCode CreateFCCeeCaloNeighbours::initialize() // retrieve decoders and other info needed for volume (ECal-HCal) connection auto decoder = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).idSpec().decoder(); - // hardcoding values (4,8) for ECal and HCal barrel systems looks quite error prone.. - // to be fixed in the future - if (m_fieldNamesSegmented[iSys] == "system" && m_fieldValuesSegmented[iSys] == 4) + if (m_connectBarrels) { - decoderECalBarrel = decoder; - // eCalThetaSize = segmentation->gridSizeTheta(); - // eCalModuleSize = 2 * M_PI / segmentation->nModules(); - // eCalModuleSize = 2 * M_PI / segmentation->phiBins(); - // eCalThetaOffset = segmentation->offsetTheta(); - // eCalPhiOffset = segmentation->offsetPhi(); - } - if (m_fieldNamesSegmented[iSys] == "system" && m_fieldValuesSegmented[iSys] == 8) - { - decoderHCalBarrel = decoder; - // hCalThetaSize = segmentation->gridSizeTheta(); - // hCalThetaOffset = segmentation->offsetTheta(); - // hCalPhiOffset = segmentation->offsetPhi(); + if (m_fieldNamesSegmented[iSys] == "system" && m_fieldValuesSegmented[iSys] == m_ecalBarrelSysId) + { + decoderECalBarrel = decoder; + ecalBarrelId = iSys; + if (segmentationType == "FCCSWGridModuleThetaMerged_k4geo") + { + ecalBarrelModuleThetaSegmentation = moduleThetaSegmentation; + } + else if (segmentationType == "FCCSWGridPhiTheta_k4geo") + { + ecalBarrelPhiThetaSegmentation = phiThetaSegmentation; + } + else + { + error() << "ECAL barrel segmentation type not handled for connectBarrels option." << endmsg; + return StatusCode::FAILURE; + } + } + + if (m_fieldNamesSegmented[iSys] == "system" && m_fieldValuesSegmented[iSys] == m_hcalBarrelSysId) + { + decoderHCalBarrel = decoder; + hcalBarrelId = iSys; + if (segmentationType == "FCCSWGridPhiTheta_k4geo") + { + hcalBarrelSegmentation = phiThetaSegmentation; + } + else + { + error() << "HCAL barrel segmentation type not handled for connectBarrels option." << endmsg; + return StatusCode::FAILURE; + } + } } + // Loop over all cells in the calorimeter and retrieve existing cellIDs and find neihbours if (segmentationType == "FCCSWGridPhiTheta_k4geo") { - // GM: code copied from RecFCChhCalorimeter, just replacing Eta->Theta - // did not review the code itself - - // Loop over all cells in the calorimeter and retrieve existing cellIDs // Loop over active layers std::vector> extrema; extrema.push_back(std::make_pair(0, m_activeVolumesNumbersSegmented[iSys] - 1)); @@ -164,15 +176,19 @@ StatusCode CreateFCCeeCaloNeighbours::initialize() (*decoder)["phi"].set(volumeId, 0); // Get number of segmentation cells within the active volume std::array numCells; - // for HCAL use theta range passed via activeVolumesTheta option - if (m_fieldValuesSegmented[iSys]==8) { - int nPhiCells = TMath::TwoPi()/phiThetaSegmentation->gridSizePhi(); + // for vlume such as HCAL for which the numberOfCells function does not work + // (because the G4 volume does not span the whole z of the detector) + // use theta range passed via activeVolumesTheta option + // if (m_fieldValuesSegmented[iSys] == m_hcalBarrelSysId) + if (m_activeVolumesTheta[iSys].size() > 0) + { + int nPhiCells = phiThetaSegmentation->phiBins(); double thetaCellSize = phiThetaSegmentation->gridSizeTheta(); double thetaOffset = phiThetaSegmentation->offsetTheta(); - double thetaMin = m_activeVolumesTheta[ilayer]; + double thetaMin = m_activeVolumesTheta[iSys][ilayer]; double thetaMax = TMath::Pi() - thetaMin; // debug - std::cout << "thetaMin, thetaMax = " << thetaMin << " " << thetaMax << std::endl; + debug() << "thetaMin, thetaMax = " << thetaMin << " " << thetaMax << endmsg; uint minThetaID = int(floor((thetaMin + 0.5 * thetaCellSize - thetaOffset) / thetaCellSize)); uint maxThetaID = int(floor((thetaMax + 0.5 * thetaCellSize - thetaOffset) / thetaCellSize)); uint nThetaCells = 1 + maxThetaID - minThetaID; @@ -180,49 +196,34 @@ StatusCode CreateFCCeeCaloNeighbours::initialize() numCells[1] = nThetaCells; numCells[2] = minThetaID; } - else { + else + { numCells = det::utils::numberOfCells(volumeId, *phiThetaSegmentation); } // extrema 1: 0, ID of last cell in phi extrema[1] = std::make_pair(0, numCells[0] - 1); - // extrema[2]: min theta ID, n theta cells + // extrema[2]: min theta ID, n theta cells extrema[2] = std::make_pair(numCells[2], numCells[1] + numCells[2] - 1); - // for layer N-1 of ECal barrel, will be used for volume connecting - /* + // for layer 0 of HCal barrel, will be used for volume connecting if ( m_fieldNamesSegmented[iSys] == "system" && - m_fieldValuesSegmented[iSys] == 4 && - ilayer == (m_activeVolumesNumbersSegmented[iSys] - 1) - ) + m_fieldValuesSegmented[iSys] == m_hcalBarrelSysId && + ilayer == 0) { - eCalLastLayer = m_activeVolumesNumbersSegmented[iSys] - 1; - // not really needed, same for all layers... unless we really restrict theta for each layer - // to the physical volume (not sure it's really needed, we can maybe allow for non-physical - // cells at edge in neighbour map, they simply won´t be added to the cluster since they do - // not exist in the list of cells) - extremaECalLastLayerPhi = std::make_pair(0, numCells[0] - 1); - extremaECalLastLayerTheta = std::make_pair(numCells[2], numCells[1] + numCells[2] - 1); + extremaHCalFirstLayerPhi = std::make_pair(0, numCells[0] - 1); + extremaHCalFirstLayerTheta = std::make_pair(numCells[2], numCells[1] + numCells[2] - 1); } - else if ( - m_fieldNamesSegmented[iSys] == "system" && - m_fieldValuesSegmented[iSys] == 8 && - m_readoutNamesSegmented[iSys] == "HCalBarrelReadout" // why is this needed? because there's also the HCAL endcap? - ) + // for layer N-1 of ECal barrel, will be used for volume connecting + if (m_fieldNamesSegmented[iSys] == "system" && + m_fieldValuesSegmented[iSys] == m_ecalBarrelSysId && + ilayer == (m_activeVolumesNumbersSegmented[iSys] - 1)) { - uint cellsTheta = ceil((2 * m_activeVolumesTheta[ilayer] - segmentation->gridSizeTheta()) / 2 / segmentation->gridSizeTheta()) * 2 + 1; // ceil( 2*m_activeVolumesRadii[ilayer] / segmentation->gridSizeEta()) ; - uint minThetaID = int(floor((-m_activeVolumesTheta[ilayer] + 0.5 * segmentation->gridSizeTheta() - segmentation->offsetTheta()) / segmentation->gridSizeTheta())); - numCells[1] = cellsTheta; - numCells[2] = minThetaID; - // for layer 0 of HCal barrel, will be used for volume connecting - if (ilayer == 0) - { - extremaHCalFirstLayerPhi = std::make_pair(0, numCells[0] - 1); - extremaHCalFirstLayerTheta = std::make_pair(numCells[2], numCells[1] + numCells[2] - 1); - extrema[2] = std::make_pair(numCells[2], numCells[1] + numCells[2] - 1); - } + eCalLastLayer = m_activeVolumesNumbersSegmented[iSys] - 1; + extremaECalLastLayerPhi = std::make_pair(0, (numCells[0] - 1)); + extremaECalLastLayerTheta = std::make_pair(numCells[2], numCells[2] + (numCells[1] - 1)); } - */ + debug() << "Layer: " << ilayer << endmsg; debug() << "Extrema[0]: " << extrema[0].first << " , " << extrema[0].second << endmsg; debug() << "Extrema[1]: " << extrema[1].first << " , " << extrema[1].second << endmsg; @@ -237,7 +238,6 @@ StatusCode CreateFCCeeCaloNeighbours::initialize() decoder->set(cellId, "phi", iphi); decoder->set(cellId, "theta", itheta + numCells[2]); // start from the minimum existing theta cell in this layer uint64_t id = cellId; - if (ilayer==7 && iphi==95) debug() << itheta+numCells[2] << " " << cellId << endmsg; map.insert(std::pair>( id, det::utils::neighbours(*decoder, {m_activeFieldNamesSegmented[iSys], "phi", "theta"}, extrema, id, {false, true, false}, m_includeDiagonalCells))); @@ -257,8 +257,6 @@ StatusCode CreateFCCeeCaloNeighbours::initialize() { dd4hep::DDSegmentation::CellID volumeId = 0; // Get VolumeID - // m_fieldValuesSegmented: in .py systemValuesModuleTheta = [4] - // m_activeFieldNamesSegmented: in .py activeFieldNamesModuleTheta = ["layer"] (*decoder)[m_fieldNamesSegmented[iSys]].set(volumeId, m_fieldValuesSegmented[iSys]); (*decoder)[m_activeFieldNamesSegmented[iSys]].set(volumeId, ilayer); (*decoder)["theta"].set(volumeId, 0); @@ -271,35 +269,19 @@ StatusCode CreateFCCeeCaloNeighbours::initialize() // extrema[2]: min theta ID, n (merged) theta cells extrema[2] = std::make_pair(numCells[2], numCells[2] + (numCells[1] - 1) * moduleThetaSegmentation->mergedThetaCells(ilayer)); // for layer N-1 of ECal barrel, will be used for volume connecting - // should 4 be systemValuesModuleTheta instead? if (ilayer == (m_activeVolumesNumbersSegmented[iSys] - 1) && m_fieldNamesSegmented[iSys] == "system" && - m_fieldValuesSegmented[iSys] == 4) - { - // eCalLastLayer = m_activeVolumesNumbersSegmented[iSys] - 1; - extremaECalLastLayerModule = std::make_pair(0, numCells[0] - 1); - extremaECalLastLayerTheta = std::make_pair(numCells[2], numCells[1] + numCells[2] - 1); - } - else if (m_fieldNamesSegmented[iSys] == "system" && - m_fieldValuesSegmented[iSys] == 8 && m_readoutNamesSegmented[iSys] == "BarHCal_Readout_phitheta") + m_fieldValuesSegmented[iSys] == m_ecalBarrelSysId) { - //// - uint cellsTheta = ceil((2 * m_activeVolumesTheta[ilayer] - segmentation->gridSizeTheta()) / 2 / segmentation->gridSizeTheta()) * 2 + 1; // ceil( 2*m_activeVolumesRadii[ilayer] / segmentation->gridSizeTheta()); - uint minThetaID = int(floor((-m_activeVolumesTheta[ilayer] + 0.5 * segmentation->gridSizeTheta() - segmentation->offsetTheta()) / segmentation->gridSizeTheta())); - numCells[1] = cellsTheta; - numCells[2] = minThetaID; - // for layer 0 of HCal barrel, will be used for volume connecting - if (ilayer == 0) - { - extremaHCalFirstLayerPhi = std::make_pair(0, numCells[0] - 1); - extremaHCalFirstLayerTheta = std::make_pair(numCells[2], numCells[1] + numCells[2] - 1); - } + eCalLastLayer = m_activeVolumesNumbersSegmented[iSys] - 1; + extremaECalLastLayerModule = std::make_pair(0, (numCells[0] - 1) * moduleThetaSegmentation->mergedModules(eCalLastLayer)); + extremaECalLastLayerTheta = std::make_pair(numCells[2], numCells[2] + (numCells[1] - 1) * moduleThetaSegmentation->mergedThetaCells(eCalLastLayer)); } debug() << "Layer: " << ilayer << endmsg; debug() << "Extrema[0]: " << extrema[0].first << " , " << extrema[0].second << endmsg; debug() << "Extrema[1]: " << extrema[1].first << " , " << extrema[1].second << endmsg; debug() << "Extrema[2]: " << extrema[2].first << " , " << extrema[2].second << endmsg; debug() << "Number of segmentation cells in (module,theta): " << numCells << endmsg; - // Loop over segmentation cells + // Loop over segmentation cells to find neighbours in ECAL for (int imodule = extrema[1].first; imodule <= extrema[1].second; imodule += moduleThetaSegmentation->mergedModules(ilayer)) { for (int itheta = extrema[2].first; itheta <= extrema[2].second; itheta += moduleThetaSegmentation->mergedThetaCells(ilayer)) @@ -349,100 +331,314 @@ StatusCode CreateFCCeeCaloNeighbours::initialize() if (m_connectBarrels) { - /* - // first check if ECAL barrel (system==5) and HCal barrel (system==8) are configured + // first check if ECAL barrel and HCal barrel are configured if (decoderECalBarrel == nullptr || decoderHCalBarrel == nullptr) { - error() << "ECAL barrel and/or HCal barrel are not configured correctly!" << endmsg; + error() << "ECAL barrel and/or HCal barrel are not configured correctly, cannot link barrels!" << endmsg; return StatusCode::FAILURE; } + // print how many cells in each dimensions will be matched info() << "ECAL layer " << eCalLastLayer << " is a neighbour of HCAL layer 0." << endmsg; - info() << "ECAL phi cells " << extremaECalLastLayerModule.first << " - " << extremaECalLastLayerModule.second - << " will be matched to HCAL cells " << extremaHCalFirstLayerPhi.first + if (ecalBarrelModuleThetaSegmentation) + { + info() << "ECAL modules " << extremaECalLastLayerModule.first << " - " << extremaECalLastLayerModule.second; + } + else + { + info() << "ECAL phi cells " << extremaECalLastLayerPhi.first << " - " << extremaECalLastLayerPhi.second; + } + info() << " will be matched to HCAL phi cells " << extremaHCalFirstLayerPhi.first << " - " << extremaHCalFirstLayerPhi.second << endmsg; info() << "ECAL theta cells " << extremaECalLastLayerTheta.first << " - " << extremaECalLastLayerTheta.second - << " will be matched to HCAL " << extremaHCalFirstLayerTheta.first + << " will be matched to HCAL theta cells " << extremaHCalFirstLayerTheta.first << " - " << extremaHCalFirstLayerTheta.second << endmsg; - std::unordered_map> thetaNeighbours; - std::unordered_map> phiNeighbours; - double hCalPhiSize = 2 * M_PI / (extremaHCalFirstLayerPhi.second - extremaHCalFirstLayerPhi.first + 1); - // loop over theta cells to match in theta - for (int iTheta = extremaHCalFirstLayerTheta.first; iTheta < extremaHCalFirstLayerTheta.second + 1; iTheta++) + dd4hep::DDSegmentation::CellID volumeIdECal = 0; + dd4hep::DDSegmentation::CellID volumeIdHCal = 0; + (*decoderECalBarrel)["system"].set(volumeIdECal, m_ecalBarrelSysId); + (*decoderECalBarrel)[m_activeFieldNamesSegmented[ecalBarrelId]].set(volumeIdECal, eCalLastLayer); + (*decoderHCalBarrel)["system"].set(volumeIdHCal, m_hcalBarrelSysId); + (*decoderHCalBarrel)[m_activeFieldNamesSegmented[hcalBarrelId]].set(volumeIdHCal, 0); + + // retrieve needed parameters + double hCalThetaSize = hcalBarrelSegmentation->gridSizeTheta(); + double hCalThetaOffset = hcalBarrelSegmentation->offsetTheta(); + double hCalPhiOffset = hcalBarrelSegmentation->offsetPhi(); + double hCalPhiSize = hcalBarrelSegmentation->gridSizePhi(); + int hCalPhiBins = hcalBarrelSegmentation->phiBins(); + double eCalThetaOffset, eCalThetaSize; + double eCalPhiOffset, eCalPhiSize; + int eCalModules(-1); + int eCalPhiBins(-1); + if (ecalBarrelModuleThetaSegmentation) + { + eCalThetaOffset = ecalBarrelModuleThetaSegmentation->offsetTheta(); + eCalThetaSize = ecalBarrelModuleThetaSegmentation->gridSizeTheta(); + eCalModules = ecalBarrelModuleThetaSegmentation->nModules(); + // for ECAL barrel with module readout, need to find out phi position of module 0 of ECal barrel last layer + // get it from volume manager + // for a normal phi grid it would suffice to use the segmentation phi(cellId) method + dd4hep::VolumeManager volman = m_geoSvc->getDetector()->volumeManager(); + auto detelement = volman.lookupDetElement(volumeIdECal); + const auto &transformMatrix = detelement.nominal().worldTransformation(); + double outGlobal[3]; + double inLocal[] = {0, 0, 0}; + transformMatrix.LocalToMaster(inLocal, outGlobal); + eCalPhiOffset = std::atan2(outGlobal[1], outGlobal[0]); + info() << "ECAL modules : " << eCalModules << " , phi of module 0 in last ECAL layer : " << eCalPhiOffset << endmsg; + eCalPhiSize = TMath::TwoPi() / eCalModules; + } + else + { + eCalThetaOffset = ecalBarrelPhiThetaSegmentation->offsetTheta(); + eCalThetaSize = ecalBarrelPhiThetaSegmentation->gridSizeTheta(); + eCalPhiOffset = ecalBarrelPhiThetaSegmentation->offsetPhi(); + eCalPhiSize = ecalBarrelPhiThetaSegmentation->gridSizePhi(); + eCalPhiBins = ecalBarrelPhiThetaSegmentation->phiBins(); + } + + dd4hep::DDSegmentation::CellID cellIdECal = volumeIdECal; + dd4hep::DDSegmentation::CellID cellIdHCal = volumeIdHCal; + + // Loop over ECAL segmentation cells to find neighbours in HCAL + if (ecalBarrelModuleThetaSegmentation) { - double lowTheta = hCalThetaOffset + iTheta * hCalThetaSize; - double highTheta = hCalThetaOffset + (iTheta + 1) * hCalThetaSize; - debug() << "HCal theta range : " << lowTheta << " - " << highTheta << endmsg; - int lowId = floor((lowTheta - 0.5 * eCalThetaSize - eCalThetaOffset) / eCalThetaSize); - int highId = floor((highTheta + 0.5 * eCalThetaSize - eCalThetaOffset) / eCalThetaSize); - debug() << "ECal theta range : " << lowId * eCalThetaSize + eCalThetaOffset << " - " - << highId * eCalThetaSize + eCalThetaOffset << endmsg; - std::vector neighbours; - for (int idThetaToAdd = lowId; idThetaToAdd <= highId; idThetaToAdd++) + for (int imodule = extremaECalLastLayerModule.first; imodule <= extremaECalLastLayerModule.second; imodule += ecalBarrelModuleThetaSegmentation->mergedModules(eCalLastLayer)) { - neighbours.push_back(det::utils::cyclicNeighbour(idThetaToAdd, extremaECalLastLayerTheta)); + (*decoderECalBarrel)["module"].set(cellIdECal, imodule); + double dphi = ecalBarrelModuleThetaSegmentation->phi(cellIdECal); + double phiVol = eCalPhiOffset + imodule * eCalPhiSize; + double phi = phiVol + dphi; + double phiMin = phi - 0.5 * eCalPhiSize * ecalBarrelModuleThetaSegmentation->mergedModules(eCalLastLayer); + double phiMax = phi + 0.5 * eCalPhiSize * ecalBarrelModuleThetaSegmentation->mergedModules(eCalLastLayer); + int minPhiHCal = (phiMin - hCalPhiOffset) / hCalPhiSize; + int maxPhiHCal = (phiMax - hCalPhiOffset) / hCalPhiSize; + for (int itheta = extremaECalLastLayerTheta.first; itheta <= extremaECalLastLayerTheta.second; itheta += ecalBarrelModuleThetaSegmentation->mergedThetaCells(eCalLastLayer)) + { + (*decoderECalBarrel)["theta"].set(cellIdECal, itheta); + double theta = ecalBarrelModuleThetaSegmentation->theta(cellIdECal); + double thetaMin = theta - eCalThetaSize * ecalBarrelModuleThetaSegmentation->mergedThetaCells(eCalLastLayer) / 2.; + double thetaMax = theta + eCalThetaSize * ecalBarrelModuleThetaSegmentation->mergedThetaCells(eCalLastLayer) / 2.; + int minThetaHCal = (thetaMin - hCalThetaOffset) / hCalThetaSize; + int maxThetaHCal = (thetaMax - hCalThetaOffset) / hCalThetaSize; + debug() << "ECAL cell: cellId = " << cellIdECal << " module = " << imodule << " theta bin = " << itheta << endmsg; + debug() << "ECAL cell: phi = " << phi << " theta = " << theta << endmsg; + debug() << "ECAL cell: thetaMin = " << thetaMin << " thetaMax = " << thetaMax << endmsg; + debug() << "ECAL cell: phiMin = " << phiMin << " phiMax = " << phiMax << endmsg; + debug() << "ECAL cell: theta bin of neighbours in HCAL layer 0 = " << minThetaHCal << " - " << maxThetaHCal << endmsg; + debug() << "ECAL cell: phi bin of neighbours in HCAL layer 0 = " << minPhiHCal << " - " << maxPhiHCal << endmsg; + bool hasNeighbours = false; + for (int iphiHCal = minPhiHCal; iphiHCal <= maxPhiHCal; iphiHCal++) + { + int phiBinHCal = iphiHCal; + // bring phi bin in 0..nbins-1 range + if (phiBinHCal < 0) + phiBinHCal += hCalPhiBins; + if (phiBinHCal >= hCalPhiBins) + phiBinHCal -= hCalPhiBins; + if (phiBinHCal < extremaHCalFirstLayerPhi.first || phiBinHCal > extremaHCalFirstLayerPhi.second) + { + warning() << "phi bin " << phiBinHCal << " out of range, skipping" << endmsg; + continue; + } + (*decoderHCalBarrel)["phi"].set(cellIdHCal, phiBinHCal); + for (int ithetaHCal = minThetaHCal; ithetaHCal <= maxThetaHCal; ithetaHCal++) + { + if (ithetaHCal < extremaHCalFirstLayerTheta.first || ithetaHCal > extremaHCalFirstLayerTheta.second) + { + warning() << "theta bin " << ithetaHCal << " out of range, skipping" << endmsg; + continue; + } + (*decoderHCalBarrel)["theta"].set(cellIdHCal, ithetaHCal); + debug() << "ECAL cell: neighbour in HCAL has cellID = " << cellIdHCal << endmsg; + map.find((uint64_t)cellIdECal)->second.push_back((uint64_t)cellIdHCal); + hasNeighbours = true; + } + } + if (hasNeighbours) + count++; + } } - debug() << "HCal theta id : " << iTheta << endmsg; - debug() << "Found ECal Neighbours in theta : " << neighbours.size() << endmsg; - for (auto id : neighbours) + } + else + { + for (int iphi = extremaECalLastLayerPhi.first; iphi <= extremaECalLastLayerPhi.second; iphi++) { - debug() << "ECal Neighbours id : " << id << endmsg; + (*decoderECalBarrel)["phi"].set(cellIdECal, iphi); + double phi = ecalBarrelPhiThetaSegmentation->phi(cellIdECal); + double phiMin = phi - 0.5 * eCalPhiSize; + double phiMax = phi + 0.5 * eCalPhiSize; + int minPhiHCal = (phiMin - hCalPhiOffset) / hCalPhiSize; + int maxPhiHCal = (phiMax - hCalPhiOffset) / hCalPhiSize; + for (int itheta = extremaECalLastLayerTheta.first; itheta <= extremaECalLastLayerTheta.second; itheta ++) + { + (*decoderECalBarrel)["theta"].set(cellIdECal, itheta); + double theta = ecalBarrelPhiThetaSegmentation->theta(cellIdECal); + double thetaMin = theta - 0.5 * eCalThetaSize; + double thetaMax = theta + 0.5 * eCalThetaSize; + int minThetaHCal = (thetaMin - hCalThetaOffset) / hCalThetaSize; + int maxThetaHCal = (thetaMax - hCalThetaOffset) / hCalThetaSize; + debug() << "ECAL cell: cellId = " << cellIdECal << " phi bin = " << iphi << " theta bin = " << itheta << endmsg; + debug() << "ECAL cell: phi = " << phi << " theta = " << theta << endmsg; + debug() << "ECAL cell: thetaMin = " << thetaMin << " thetaMax = " << thetaMax << endmsg; + debug() << "ECAL cell: phiMin = " << phiMin << " phiMax = " << phiMax << endmsg; + debug() << "ECAL cell: theta bin of neighbours in HCAL layer 0 = " << minThetaHCal << " - " << maxThetaHCal << endmsg; + debug() << "ECAL cell: phi bin of neighbours in HCAL layer 0 = " << minPhiHCal << " - " << maxPhiHCal << endmsg; + bool hasNeighbours = false; + for (int iphiHCal = minPhiHCal; iphiHCal <= maxPhiHCal; iphiHCal++) + { + int phiBinHCal = iphiHCal; + // bring phi bin in 0..nbins-1 range + if (phiBinHCal < 0) + phiBinHCal += hCalPhiBins; + if (phiBinHCal >= hCalPhiBins) + phiBinHCal -= hCalPhiBins; + if (phiBinHCal < extremaHCalFirstLayerPhi.first || phiBinHCal > extremaHCalFirstLayerPhi.second) + { + warning() << "phi bin " << phiBinHCal << " out of range, skipping" << endmsg; + continue; + } + (*decoderHCalBarrel)["phi"].set(cellIdHCal, phiBinHCal); + for (int ithetaHCal = minThetaHCal; ithetaHCal <= maxThetaHCal; ithetaHCal++) + { + if (ithetaHCal < extremaHCalFirstLayerTheta.first || ithetaHCal > extremaHCalFirstLayerTheta.second) + { + warning() << "theta bin " << ithetaHCal << " out of range, skipping" << endmsg; + continue; + } + (*decoderHCalBarrel)["theta"].set(cellIdHCal, ithetaHCal); + debug() << "ECAL cell: neighbour in HCAL has cellID = " << cellIdHCal << endmsg; + map.find((uint64_t)cellIdECal)->second.push_back((uint64_t)cellIdHCal); + hasNeighbours = true; + } + } + if (hasNeighbours) + count++; + } } - thetaNeighbours.insert(std::pair>(iTheta, neighbours)); } - // loop over phi and find which phi cells to add - for (int iPhi = 0; iPhi < extremaHCalFirstLayerPhi.second + 1; iPhi++) + // Loop over HCAL segmentation cells to find neighbours in ECAL + // Loop on phi + for (int iphi = extremaHCalFirstLayerPhi.first; iphi <= extremaHCalFirstLayerPhi.second; iphi++) { - double lowPhi = hCalPhiOffset + iPhi * hCalPhiSize; - double highPhi = hCalPhiOffset + (iPhi + 1) * hCalPhiSize; - debug() << "HCal phi range : " << lowPhi << " - " << highPhi << endmsg; - int lowId = floor((lowPhi - 0.5 * eCalModuleSize - eCalPhiOffset) / eCalModuleSize); - int highId = floor((highPhi + 0.5 * eCalModuleSize - eCalPhiOffset) / eCalModuleSize); - debug() << "ECal phi range : " << lowId * eCalModuleSize + eCalPhiOffset << " - " - << highId * eCalModuleSize + eCalPhiOffset << endmsg; - std::vector neighbours; - for (int idPhiToAdd = lowId; idPhiToAdd <= highId; idPhiToAdd++) + // determine phi extent of HCAL cell + (*decoderHCalBarrel)["phi"].set(cellIdHCal, iphi); + double phi = hcalBarrelSegmentation->phi(cellIdHCal); + double phiMin = phi - 0.5 * hCalPhiSize; + double phiMax = phi + 0.5 * hCalPhiSize; + + // find ECAL barrel modules corresponding to this phi range + int minModuleECal(-1), maxModuleECal(-1), minPhiECal(-1), maxPhiECal(-1); + if (ecalBarrelModuleThetaSegmentation) { - neighbours.push_back(det::utils::cyclicNeighbour(idPhiToAdd, extremaECalLastLayerModule)); + minModuleECal = (int)((phiMin - eCalPhiOffset) / eCalPhiSize); + if (minModuleECal < 0) + minModuleECal += eCalModules; // need module to be >=0 so that subtracting module%mergedModules gives the correct result + minModuleECal -= minModuleECal % ecalBarrelModuleThetaSegmentation->mergedModules(eCalLastLayer); + maxModuleECal = (int)((phiMax - eCalPhiOffset) / eCalPhiSize); + if (maxModuleECal < 0) + maxModuleECal += eCalModules; + maxModuleECal -= maxModuleECal % ecalBarrelModuleThetaSegmentation->mergedModules(eCalLastLayer); + // due to ciclic behaviour of modules, one could have e.g min = 1534 and max = 2 (for N=1536) + if (maxModuleECal < minModuleECal) + maxModuleECal += eCalModules; } - debug() << "HCal phi id : " << iPhi << endmsg; - debug() << "Found ECal Neighbours in phi : " << neighbours.size() << endmsg; - for (auto id : neighbours) + else { - debug() << "ECal Neighbours id : " << id << endmsg; + minPhiECal = (int)((phiMin - eCalPhiOffset) / eCalPhiSize); + maxPhiECal = (int)((phiMax - eCalPhiOffset) / eCalPhiSize); } - phiNeighbours.insert(std::pair>(iPhi, neighbours)); - } - // add neighbours to both ecal cell and hcal cells - dd4hep::DDSegmentation::CellID ecalCellId = 0; - dd4hep::DDSegmentation::CellID hcalCellId = 0; - (*decoderECalBarrel)["system"].set(ecalCellId, 5); - (*decoderECalBarrel)[m_activeFieldNamesSegmented[0]].set(ecalCellId, eCalLastLayer); - (*decoderHCalBarrel)["system"].set(hcalCellId, 8); - (*decoderHCalBarrel)[m_activeFieldNamesSegmented[1]].set(hcalCellId, 0); - // loop over segmented hcal cells - for (auto iThetaHCal : thetaNeighbours) - { - (*decoderHCalBarrel)["theta"].set(hcalCellId, iThetaHCal.first); - for (auto iPhiHCal : phiNeighbours) + + // Loop on theta + for (int itheta = extremaHCalFirstLayerTheta.first; itheta <= extremaHCalFirstLayerTheta.second; itheta++) { - (*decoderHCalBarrel)["phi"].set(hcalCellId, iPhiHCal.first); - for (auto iTheta : iThetaHCal.second) + // determine theta extent of HCAL cell + (*decoderHCalBarrel)["theta"].set(cellIdHCal, itheta); + double theta = hcalBarrelSegmentation->theta(cellIdHCal); + double thetaMin = theta - 0.5 * hCalThetaSize; + double thetaMax = theta + 0.5 * hCalThetaSize; + // find ECAL barrel theta bins corresponding to this theta range + int minThetaECal = (thetaMin - eCalThetaOffset) / eCalThetaSize; + int maxThetaECal = (thetaMax - eCalThetaOffset) / eCalThetaSize; + if (ecalBarrelModuleThetaSegmentation) + { + minThetaECal -= minThetaECal % ecalBarrelModuleThetaSegmentation->mergedThetaCells(eCalLastLayer); + maxThetaECal -= maxThetaECal % ecalBarrelModuleThetaSegmentation->mergedThetaCells(eCalLastLayer); + } + // print out some debug info + debug() << "HCAL cell: cellId = " << cellIdHCal << " phi bin = " << iphi << " theta bin = " << itheta << endmsg; + debug() << "HCAL cell: phi = " << phi << " theta = " << theta << endmsg; + debug() << "HCAL cell: thetaMin = " << thetaMin << " thetaMax = " << thetaMax << endmsg; + debug() << "HCAL cell: phiMin = " << phiMin << " phiMax = " << phiMax << endmsg; + debug() << "HCAL cell: theta bin of neighbours in ECAL layer max = " << minThetaECal << " - " << maxThetaECal << endmsg; + if (ecalBarrelModuleThetaSegmentation) + { + debug() << "HCAL cell: module of neighbours in ECAL layer max = " << minModuleECal << " - " << maxModuleECal << endmsg; + } + else { - (*decoderECalBarrel)["theta"].set(ecalCellId, iTheta); - for (auto iPhi : iPhiHCal.second) + debug() << "HCAL cell: phi bin of neighbours in ECAL layer max = " << minPhiECal << " - " << maxPhiECal << endmsg; + } + + // add the neighbours if within the acceptance of the layer + bool hasNeighbours = false; + int thetaStep = (ecalBarrelModuleThetaSegmentation) ? ecalBarrelModuleThetaSegmentation->mergedThetaCells(eCalLastLayer) : 1; + for (int ithetaECal = minThetaECal; ithetaECal <= maxThetaECal; ithetaECal += thetaStep) + { + if (ithetaECal < extremaECalLastLayerTheta.first || ithetaECal > extremaECalLastLayerTheta.second) { - (*decoderECalBarrel)["phi"].set(ecalCellId, iPhi); - map.find(hcalCellId)->second.push_back(ecalCellId); - map.find(ecalCellId)->second.push_back(hcalCellId); - count++; + debug() << "theta bin " << ithetaECal << " out of range, skipping" << endmsg; + continue; } + (*decoderECalBarrel)["theta"].set(cellIdECal, ithetaECal); + if (ecalBarrelModuleThetaSegmentation) + { + for (int iModule = minModuleECal; iModule <= maxModuleECal; iModule += ecalBarrelModuleThetaSegmentation->mergedModules(eCalLastLayer)) + { + int module = iModule; + if (module < 0) + module += eCalModules; + if (module >= eCalModules) + module -= eCalModules; + if (module < extremaECalLastLayerModule.first || module > extremaECalLastLayerModule.second) + { + warning() << "module " << module << " out of range, skipping" << endmsg; + continue; + } + (*decoderECalBarrel)["module"].set(cellIdECal, module); + debug() << "HCAL cell: neighbour in ECAL has cellID = " << cellIdECal << endmsg; + map.find((uint64_t)cellIdHCal)->second.push_back((uint64_t)cellIdECal); + hasNeighbours = true; + } + } + else + { + for (int iphiECal = minPhiECal; iphiECal <= maxPhiECal; iphiECal++) + { + int phiBinECal = iphiECal; + // bring phi bin in 0..nbins-1 range + if (phiBinECal < 0) + phiBinECal += eCalPhiBins; + if (phiBinECal >= eCalPhiBins) + phiBinECal -= eCalPhiBins; + if (phiBinECal < extremaECalLastLayerPhi.first || phiBinECal > extremaECalLastLayerPhi.second) + { + warning() << "phi bin " << phiBinECal << " out of range, skipping" << endmsg; + continue; + } + (*decoderECalBarrel)["phi"].set(cellIdECal, phiBinECal); + debug() << "HCAL cell: neighbour in ECAL has cellID = " << cellIdECal << endmsg; + map.find((uint64_t)cellIdHCal)->second.push_back((uint64_t)cellIdECal); + hasNeighbours = true; + } + } + if (hasNeighbours) + count++; } } } - */ } // end of connectBarrels diff --git a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.h b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.h index 2f60de9b..9c256e0a 100644 --- a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.h +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNeighbours.h @@ -50,23 +50,24 @@ class CreateFCCeeCaloNeighbours : public extends1 { /// Pointer to the geometry service SmartIF m_geoSvc; - /// Names of the detector readout for volumes with theta-module segmentation - Gaudi::Property> m_readoutNamesSegmented{this, "readoutNames", {"ECalBarrelModuleThetaMerged"}}; + /// Names of the detector readout for the volumes + Gaudi::Property> m_readoutNamesSegmented{this, "readoutNames", {"ECalBarrelModuleThetaMerged", "BarHCal_Readout_phitheta"}}; /// Name of the fields describing the segmented volume - Gaudi::Property> m_fieldNamesSegmented{this, "systemNames", {"system"}}; + Gaudi::Property> m_fieldNamesSegmented{this, "systemNames", {"system", "system"}}; /// Values of the fields describing the segmented volume - Gaudi::Property> m_fieldValuesSegmented{this, "systemValues", {5}}; - /// Names of the active volume in geometry along radial axis (e.g. layer), the others are "module"/"phi", "theta" - Gaudi::Property> m_activeFieldNamesSegmented{this, "activeFieldNames", {"layer"}}; - /// Number of layers in the segmented volume - Gaudi::Property> m_activeVolumesNumbersSegmented{this, "activeVolumesNumbers", {12}}; - // Radii of layers in the segmented volume - Gaudi::Property> m_activeVolumesTheta{this, "activeVolumesTheta"}; - /// Whether to create the geant4 geometry or not + Gaudi::Property> m_fieldValuesSegmented{this, "systemValues", {4, 8}}; + /// Names of the active volume in geometry along radial axis (e.g. layer), the others are "module" or "phi", "theta" + Gaudi::Property> m_activeFieldNamesSegmented{this, "activeFieldNames", {"layer", "layer"}}; + /// Number of layers in the segmented volumes + Gaudi::Property> m_activeVolumesNumbersSegmented{this, "activeVolumesNumbers", {12, 13}}; + // Theta ranges of layers in the segmented volumes + Gaudi::Property>> m_activeVolumesTheta{this, "activeVolumesTheta"}; + /// Whether to consider diagonal cells as neighbours or not Gaudi::Property m_includeDiagonalCells{this, "includeDiagonalCells", false, "If True will consider also diagonal neighbours in volumes with theta-module segmentation"}; - - /// Name of output file - std::string m_outputFileName; + + // System ID of ECAL and HCAL barrels + Gaudi::Property m_ecalBarrelSysId{this, "ecalBarrelSysId", 4}; + Gaudi::Property m_hcalBarrelSysId{this, "hcalBarrelSysId", 8}; // For combination of barrels: flag if ECal and HCal barrels should be merged Gaudi::Property m_connectBarrels{this, "connectBarrels", true}; @@ -78,6 +79,9 @@ class CreateFCCeeCaloNeighbours : public extends1 { Gaudi::Property m_hCalRinner{this, "hCalRinner", 2850}; // For combination of barrels: offset of HCal modules in phi (lower edge) Gaudi::Property m_hCalPhiOffset{this, "hCalPhiOffset"}; + + /// Name of output file + std::string m_outputFileName; }; #endif /* RECFCCEECALORIMETER_CREATEFCCEECALONEIGHBOURS_H */ diff --git a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp index ac751338..8a4562f2 100644 --- a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp @@ -10,109 +10,230 @@ DECLARE_COMPONENT(CreateFCCeeCaloNoiseLevelMap) -CreateFCCeeCaloNoiseLevelMap::CreateFCCeeCaloNoiseLevelMap(const std::string& aName, ISvcLocator* aSL) - : base_class(aName, aSL) { - declareProperty("ECalBarrelNoiseTool", m_ecalBarrelNoiseTool, "Handle for the cells noise tool of Barrel ECal"); - declareProperty("HCalBarrelNoiseTool", m_hcalBarrelNoiseTool, "Handle for the cells noise tool of Barrel HCal"); - declareProperty( "outputFileName", m_outputFileName, "Name of the output file"); +CreateFCCeeCaloNoiseLevelMap::CreateFCCeeCaloNoiseLevelMap(const std::string &aName, ISvcLocator *aSL) + : base_class(aName, aSL) +{ + declareProperty("ECalBarrelNoiseTool", m_ecalBarrelNoiseTool, "Handle for the cell noise tool of Barrel ECal"); + declareProperty("HCalBarrelNoiseTool", m_hcalBarrelNoiseTool, "Handle for the cell noise tool of Barrel HCal"); + declareProperty("outputFileName", m_outputFileName, "Name of the output file"); } CreateFCCeeCaloNoiseLevelMap::~CreateFCCeeCaloNoiseLevelMap() {} -StatusCode CreateFCCeeCaloNoiseLevelMap::initialize() { +StatusCode CreateFCCeeCaloNoiseLevelMap::initialize() +{ // Initialize necessary Gaudi components - if (Service::initialize().isFailure()) { - error() << "Unable to initialize Service()" << endmsg; + if (Service::initialize().isFailure()) + { return StatusCode::FAILURE; } m_geoSvc = service("GeoSvc"); - if (!m_geoSvc) { + 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_ecalBarrelNoiseTool.retrieve()) { + if (!m_ecalBarrelNoiseTool.retrieve()) + { error() << "Unable to retrieve the ECAL noise tool!!!" << endmsg; return StatusCode::FAILURE; } - if (!m_hcalBarrelNoiseTool.empty()){ - if (!m_hcalBarrelNoiseTool.retrieve()) { + if (!m_hcalBarrelNoiseTool.empty()) + { + if (!m_hcalBarrelNoiseTool.retrieve()) + { error() << "Unable to retrieve the HCAL noise tool!!!" << endmsg; return StatusCode::FAILURE; } } - - std::unordered_map> map; - for (uint iSys = 0; iSys < m_readoutNamesSegmented.size(); iSys++) { + std::unordered_map> map; + + for (uint iSys = 0; iSys < m_readoutNamesSegmented.size(); iSys++) + { // Check if readouts exist info() << "Readout: " << m_readoutNamesSegmented[iSys] << endmsg; - if (m_geoSvc->getDetector()->readouts().find(m_readoutNamesSegmented[iSys]) == m_geoSvc->getDetector()->readouts().end()) { + if (m_geoSvc->getDetector()->readouts().find(m_readoutNamesSegmented[iSys]) == m_geoSvc->getDetector()->readouts().end()) + { error() << "Readout <<" << m_readoutNamesSegmented[iSys] << ">> does not exist." << endmsg; return StatusCode::FAILURE; } + // get segmentation - dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo* segmentation; - segmentation = dynamic_cast( - m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).segmentation().segmentation()); - if (segmentation == nullptr) { - error() << "There is no Module-theta segmentation!!!!" << endmsg; + dd4hep::DDSegmentation::Segmentation *aSegmentation = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).segmentation().segmentation(); + std::string segmentationType = aSegmentation->type(); + info() << "Segmentation type : " << segmentationType << endmsg; + + dd4hep::DDSegmentation::GridTheta_k4geo *segmentation = nullptr; + dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo *phiThetaSegmentation = nullptr; + dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo *moduleThetaSegmentation = nullptr; + if (segmentationType == "FCCSWGridModuleThetaMerged_k4geo") + { + segmentation = dynamic_cast(aSegmentation); + moduleThetaSegmentation = dynamic_cast(aSegmentation); + } + else if (segmentationType == "FCCSWGridPhiTheta_k4geo") + { + segmentation = dynamic_cast(aSegmentation); + phiThetaSegmentation = dynamic_cast(aSegmentation); + } + else + { + error() << "Segmentation type not handled." << endmsg; + return StatusCode::FAILURE; + } + + if (segmentation == nullptr || (phiThetaSegmentation == nullptr && moduleThetaSegmentation == nullptr)) + { + error() << "Unable to cast segmentation pointer!!!!" << endmsg; return StatusCode::FAILURE; } - info() << "FCCSWGridModuleThetaMerged: size in Theta " << segmentation->gridSizeTheta() << " , bins in Module " << segmentation->nModules() - << endmsg; - info() << "FCCSWGridModuleThetaMerged: offset in Theta " << segmentation->offsetTheta() << endmsg; + info() << "Segmentation: size in Theta " << segmentation->gridSizeTheta() << endmsg; + info() << "Segmentation: offset in Theta " << segmentation->offsetTheta() << endmsg; + if (segmentationType == "FCCSWGridModuleThetaMerged_k4geo") + { + info() << "Segmentation: bins in Module " << moduleThetaSegmentation->nModules() << endmsg; + } + else if (segmentationType == "FCCSWGridPhiTheta_k4geo") + { + info() << "Segmentation: size in Phi " << phiThetaSegmentation->gridSizePhi() << endmsg; + info() << "Segmentation: offset in Phi " << phiThetaSegmentation->offsetPhi() << endmsg; + } - auto decoder = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).idSpec().decoder(); // Loop over all cells in the calorimeter and retrieve existing cellIDs // Loop over active layers + auto decoder = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).idSpec().decoder(); std::vector> extrema; + // extrema[0]: (0, nLayers-1) extrema.push_back(std::make_pair(0, m_activeVolumesNumbersSegmented[iSys] - 1)); + // extrema[1] and extrema[2] will contain theta and phi/module ranges, will be set later + // as they (theta) depend on layer extrema.push_back(std::make_pair(0, 0)); extrema.push_back(std::make_pair(0, 0)); - for (unsigned int ilayer = 0; ilayer < m_activeVolumesNumbersSegmented[iSys]; ilayer++) { - dd4hep::DDSegmentation::CellID volumeId = 0; - // Get VolumeID - (*decoder)[m_fieldNamesSegmented[iSys]].set(volumeId, m_fieldValuesSegmented[iSys]); - (*decoder)[m_activeFieldNamesSegmented[iSys]].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, *segmentation); - // Range of module ID - extrema[1] = std::make_pair(0, (numCells[0] - 1)*segmentation->mergedModules(ilayer)); - // Range and min of theta ID - // for HCAL use alternative calculation - if (m_fieldNamesSegmented[iSys] == "system" && - m_fieldValuesSegmented[iSys] == m_hcalBarrelSysId){ - uint cellsTheta = ceil(( 2*m_activeVolumesTheta[ilayer] - segmentation->gridSizeTheta() ) / 2 / segmentation->gridSizeTheta()) * 2 + 1; //ceil( 2*m_activeVolumesRadii[ilayer] / segmentation->gridSizeTheta()); - uint minThetaID = int(floor(( - m_activeVolumesTheta[ilayer] + 0.5 * segmentation->gridSizeTheta() - segmentation->offsetTheta()) / segmentation->gridSizeTheta())); - numCells[1]=cellsTheta; - numCells[2]=minThetaID; + for (unsigned int ilayer = 0; ilayer < m_activeVolumesNumbersSegmented[iSys]; ilayer++) + { + if (segmentationType == "FCCSWGridPhiTheta_k4geo") + { + // for ECAL, G4 volumes extend along the full z so + // they can be used - once obtained from the cellID - to + // determine the full theta range of the layer. + // For HCal the G4 volume containing the cell spans only + // a subrange of the full z range of the HCal layer - one + // would need to find the HCalLayerVol volume to determine + // the theta range from it. + // Thus for HCAL we use the theta range passed via activeVolumesTheta + + // Set system and layer of volume ID + dd4hep::DDSegmentation::CellID volumeId = 0; + (*decoder)[m_fieldNamesSegmented[iSys]].set(volumeId, m_fieldValuesSegmented[iSys]); + (*decoder)[m_activeFieldNamesSegmented[iSys]].set(volumeId, ilayer); + (*decoder)["theta"].set(volumeId, 0); + (*decoder)["phi"].set(volumeId, 0); + // Get number of segmentation cells within the active volume + std::array numCells; + // for volumes (such as Allegro's HCAL) for which theta range + // is passed via activeVolumesTheta option then use that instead of + // numberOfCells helper function + if (m_activeVolumesTheta[iSys].size()>0) + { + int nPhiCells = TMath::TwoPi() / phiThetaSegmentation->gridSizePhi(); + double thetaCellSize = phiThetaSegmentation->gridSizeTheta(); + double thetaOffset = phiThetaSegmentation->offsetTheta(); + double thetaMin = m_activeVolumesTheta[iSys][ilayer]; + double thetaMax = TMath::Pi() - thetaMin; + uint minThetaID = int(floor((thetaMin + 0.5 * thetaCellSize - thetaOffset) / thetaCellSize)); + uint maxThetaID = int(floor((thetaMax + 0.5 * thetaCellSize - thetaOffset) / thetaCellSize)); + uint nThetaCells = 1 + maxThetaID - minThetaID; + numCells[0] = nPhiCells; + numCells[1] = nThetaCells; + numCells[2] = minThetaID; + } + else + { + numCells = det::utils::numberOfCells(volumeId, *phiThetaSegmentation); + } + // extrema[1]: 0, ID of last cell in phi + extrema[1] = std::make_pair(0, numCells[0] - 1); + // extrema[2]: min theta ID, ID of last cell in theta + extrema[2] = std::make_pair(numCells[2], numCells[1] + numCells[2] - 1); + debug() << "Layer: " << ilayer << endmsg; + debug() << "Extrema[0]: " << extrema[0].first << " , " << extrema[0].second << endmsg; + debug() << "Extrema[1]: " << extrema[1].first << " , " << extrema[1].second << endmsg; + debug() << "Extrema[2]: " << extrema[2].first << " , " << extrema[2].second << endmsg; + debug() << "Number of segmentation cells in phi, in theta, and min theta ID, : " << numCells << endmsg; + // Loop over segmentation cells + for (unsigned int iphi = 0; iphi < numCells[0]; iphi++) + { + for (unsigned int itheta = 0; itheta < numCells[1]; itheta++) + { + dd4hep::DDSegmentation::CellID cellId = volumeId; + decoder->set(cellId, "phi", iphi); + decoder->set(cellId, "theta", itheta + numCells[2]); // start from the minimum existing theta cell in this layer + uint64_t id = cellId; + double noise = 0.; + double noiseOffset = 0.; + if (m_fieldValuesSegmented[iSys] == m_ecalBarrelSysId) + { + noise = m_ecalBarrelNoiseTool->getNoiseConstantPerCell(id); + noiseOffset = m_ecalBarrelNoiseTool->getNoiseOffsetPerCell(id); + } + else if (m_fieldValuesSegmented[iSys] == m_hcalBarrelSysId) + { + noise = m_hcalBarrelNoiseTool->getNoiseConstantPerCell(id); + noiseOffset = m_hcalBarrelNoiseTool->getNoiseOffsetPerCell(id); + } + else + { + warning() << "Unexpected system value for phi-theta readout " << m_fieldValuesSegmented[iSys] + << ", setting noise to 0.0" << endmsg; + } + map.insert(std::pair>(id, std::make_pair(noise, noiseOffset))); + } + } } - extrema[2] = std::make_pair(numCells[2], numCells[2] + (numCells[1] - 1)*segmentation->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*segmentation->mergedModules(ilayer)); - decoder->set(cellId, "theta", numCells[2] + itheta*segmentation->mergedThetaCells(ilayer)); // start from the minimum existing theta cell in this layer - uint64_t id = cellId; - double noise = 0.; - double noiseOffset = 0.; - if (m_fieldValuesSegmented[iSys] == m_hcalBarrelSysId){ - noise = m_hcalBarrelNoiseTool->getNoiseConstantPerCell(id); - noiseOffset = m_hcalBarrelNoiseTool->getNoiseOffsetPerCell(id); - } else if (m_fieldValuesSegmented[iSys] == m_ecalBarrelSysId){ - noise = m_ecalBarrelNoiseTool->getNoiseConstantPerCell(id); - noiseOffset = m_ecalBarrelNoiseTool->getNoiseOffsetPerCell(id); - } - map.insert( std::pair >(id, std::make_pair(noise, noiseOffset) ) ); + else if (segmentationType == "FCCSWGridModuleThetaMerged_k4geo") + { + // Set system and layer of volume ID + dd4hep::DDSegmentation::CellID volumeId = 0; + (*decoder)[m_fieldNamesSegmented[iSys]].set(volumeId, m_fieldValuesSegmented[iSys]); + (*decoder)[m_activeFieldNamesSegmented[iSys]].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 + uint64_t id = cellId; + double noise = 0.; + double noiseOffset = 0.; + if (m_fieldValuesSegmented[iSys] == m_ecalBarrelSysId) + { + noise = m_ecalBarrelNoiseTool->getNoiseConstantPerCell(id); + noiseOffset = m_ecalBarrelNoiseTool->getNoiseOffsetPerCell(id); + } + else + { + warning() << "Unexpected system value for module-theta readout " << m_fieldValuesSegmented[iSys] + << ", setting noise to 0.0" << endmsg; + } + map.insert(std::pair>(id, std::make_pair(noise, noiseOffset))); + } } } } @@ -127,7 +248,8 @@ StatusCode CreateFCCeeCaloNoiseLevelMap::initialize() { tree.Branch("cellId", &saveCellId, "cellId/l"); tree.Branch("noiseLevel", &saveNoiseLevel); tree.Branch("noiseOffset", &saveNoiseOffset); - for (const auto& item : map) { + for (const auto &item : map) + { saveCellId = item.first; saveNoiseLevel = item.second.first; saveNoiseOffset = item.second.second; diff --git a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.h b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.h index 44f7eb4e..b682094a 100644 --- a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.h +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.h @@ -40,9 +40,11 @@ class CreateFCCeeCaloNoiseLevelMap : public extends1 { /// ideally could replace with vector of handles (so that we can also have ecal endcap and so on) /// Handle for the cells noise tool in ECal ToolHandle m_ecalBarrelNoiseTool{"ReadNoiseFromFileTool", this}; - Gaudi::Property m_ecalBarrelSysId{this, "ecalBarrelSysId", 5}; /// Handle for the cells noise tool in HCal ToolHandle m_hcalBarrelNoiseTool{"ReadNoiseFromFileTool", this}; + + // System ID of ECAL and HCAL barrels + Gaudi::Property m_ecalBarrelSysId{this, "ecalBarrelSysId", 4}; Gaudi::Property m_hcalBarrelSysId{this, "hcalBarrelSysId", 8}; /// Names of the detector readout for the volumes @@ -52,11 +54,11 @@ class CreateFCCeeCaloNoiseLevelMap : public extends1 { /// Values of the fields describing the segmented volume Gaudi::Property> m_fieldValuesSegmented{this, "systemValues", {4, 8}}; /// Names of the active volume in geometry along radial axis (e.g. layer), the others are "module" or "phi", "theta" - Gaudi::Property> m_activeFieldNamesSegmented{this, "activeFieldNames", {"layer"}}; - /// Number of layers in the segmented volume + Gaudi::Property> m_activeFieldNamesSegmented{this, "activeFieldNames", {"layer", "layer"}}; + /// Number of layers in the segmented volumes Gaudi::Property> m_activeVolumesNumbersSegmented{this, "activeVolumesNumbers", {12, 13}}; - // Radii of layers in the segmented volume - Gaudi::Property> m_activeVolumesTheta{this, "activeVolumesTheta"}; + // Theta ranges of layers in the segmented volumes (if needed) + Gaudi::Property>> m_activeVolumesTheta{this, "activeVolumesTheta"}; /// Name of output file std::string m_outputFileName;