diff --git a/RecCalorimeter/src/components/Merge2CaloHitCollections.cpp b/RecCalorimeter/src/components/Merge2CaloHitCollections.cpp new file mode 100644 index 00000000..c2a7d707 --- /dev/null +++ b/RecCalorimeter/src/components/Merge2CaloHitCollections.cpp @@ -0,0 +1,59 @@ +#include "Merge2CaloHitCollections.h" +#include + + +DECLARE_COMPONENT(Merge2CaloHitCollections) + + +Merge2CaloHitCollections::Merge2CaloHitCollections( + const std::string& name, + ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { + declareProperty("inHitsA", m_inCollA, "First calo hit collection (input)"); + declareProperty("inHitsB", m_inCollB, "Second calo hit collection (input)"); + declareProperty("outHits", m_outColl, "Resulting hit collection (output)"); +} + + +StatusCode Merge2CaloHitCollections::initialize() { + { + StatusCode sc = GaudiAlgorithm::initialize(); + if (sc.isFailure()) { + return sc; + } + } + + debug() << "Merge 2 calo hit collections initialized!" << endmsg; + + return StatusCode::SUCCESS; +} + + +StatusCode Merge2CaloHitCollections::execute() { + // Create a new empty calo hits collection + edm4hep::CalorimeterHitCollection* outColl = m_outColl.createAndPut(); + + // Clone hits from the first collection + const edm4hep::CalorimeterHitCollection* inCollA = m_inCollA.get(); + debug() << "Cloning " << inCollA->size() << " from the \"" + << m_inCollA.fullKey().key() << "\" collection." << endmsg; + for (const auto& hit : *inCollA) { + outColl->push_back(hit.clone()); + } + + // Clone hits from the Second collection + const edm4hep::CalorimeterHitCollection* inCollB = m_inCollB.get(); + debug() << "Cloning " << inCollB->size() << " from the \"" + << m_inCollB.fullKey().key() << "\" collection." << endmsg; + for (const auto& hit : *inCollB) { + outColl->push_back(hit.clone()); + } + + debug() << "Output collection \"" << m_outColl.fullKey().key() + << "\" contains " << outColl->size() << " calo hits." << endmsg; + return StatusCode::SUCCESS; +} + + +StatusCode Merge2CaloHitCollections::finalize() { + return GaudiAlgorithm::finalize(); +} diff --git a/RecCalorimeter/src/components/Merge2CaloHitCollections.h b/RecCalorimeter/src/components/Merge2CaloHitCollections.h new file mode 100644 index 00000000..09606399 --- /dev/null +++ b/RecCalorimeter/src/components/Merge2CaloHitCollections.h @@ -0,0 +1,53 @@ +#ifndef RECCALORIMETER_MERGE_2_CALO_HIT_COLLECTIONS_H +#define RECCALORIMETER_MERGE_2_CALO_HIT_COLLECTIONS_H + +// k4FWCore +#include "k4FWCore/DataHandle.h" + +// Gaudi +#include "GaudiAlg/GaudiAlgorithm.h" + +// Datamodel +#include "edm4hep/CalorimeterHitCollection.h" + +class IGeoSvc; + +/** @class Merge2CaloHitCollections + * + * Algorithm for merging two collections together + * Use-case: + * Merging multiple calorimeter collections together. + * (ECAL barrel, HCAL barrel+extended barrel, ECAL + HCAL endcaps, + * ECAL + HCAL forward calorimeter). + * + * @author Juraj Smiesko + * @date 2023-11-06 + * + */ + +class Merge2CaloHitCollections : public GaudiAlgorithm { + +public: + Merge2CaloHitCollections(const std::string& name, ISvcLocator* svcLoc); + + StatusCode initialize(); + + StatusCode execute(); + + StatusCode finalize(); + +private: + /// Handle for the first hit collection (input) + DataHandle m_inCollA{ + "inHitsA", Gaudi::DataHandle::Reader, this}; + + /// Handle for the second hit collection (input) + DataHandle m_inCollB{ + "inHitsB", Gaudi::DataHandle::Reader, this}; + + /// Handle for the resulting hit collection (output) + DataHandle m_outColl{ + "outHits", Gaudi::DataHandle::Writer, this}; +}; + +#endif /* RECCALORIMETER_MERGE_2_CALO_HIT_COLLECTIONS_H */ diff --git a/RecCalorimeter/src/components/TopoCaloNeighbours.h b/RecCalorimeter/src/components/TopoCaloNeighbours.h index 9cd70f61..7179efac 100644 --- a/RecCalorimeter/src/components/TopoCaloNeighbours.h +++ b/RecCalorimeter/src/components/TopoCaloNeighbours.h @@ -23,12 +23,12 @@ class TopoCaloNeighbours : public GaudiTool, virtual public ICaloReadNeighboursM public: TopoCaloNeighbours(const std::string& type, const std::string& name, const IInterface* parent); virtual ~TopoCaloNeighbours() = default; - + /** Read a map of cellIDs to vector of cellIDs (neighbours). */ virtual StatusCode initialize() final; virtual StatusCode finalize() final; - + /** Function to be called for the neighbours of a cell. * @param[in] aCellId, cellid of the cell of interest. * @return vector of cellIDs, corresponding to the cells neighbours. diff --git a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp index ebeb9d8e..aef71bcb 100644 --- a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp +++ b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp @@ -1,11 +1,26 @@ #include "CaloTopoClusterFCCee.h" #include "../../../RecCalorimeter/src/components/NoiseCaloCellsFromFileTool.h" -// FCCSW -#include "DetCommon/DetUtils.h" +// std +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Gaudi +#include + +// k4FWCore #include "k4Interface/IGeoSvc.h" -// datamodel +// Datamodel #include "edm4hep/Cluster.h" #include "edm4hep/ClusterCollection.h" #include "edm4hep/CalorimeterHit.h" @@ -15,438 +30,452 @@ #include "DD4hep/Detector.h" #include "DD4hep/Readout.h" -#include -#include -#include -#include -#include -#include -#include DECLARE_COMPONENT(CaloTopoClusterFCCee) + CaloTopoClusterFCCee::CaloTopoClusterFCCee(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) { - declareProperty("TopoClusterInput", m_inputTool, "Handle for input map of cells"); - declareProperty("noiseTool", m_noiseTool, "Handle for the cells noise tool"); - declareProperty("neigboursTool", m_neighboursTool, "Handle for tool to retrieve cell neighbours"); - declareProperty("positionsECalBarrelTool", m_cellPositionsECalBarrelTool, - "Handle for tool to retrieve cell positions in ECal Barrel"); - declareProperty("positionsHCalBarrelTool", m_cellPositionsHCalBarrelTool, - "Handle for tool to retrieve cell positions in HCal Barrel"); - declareProperty("positionsHCalBarrelNoSegTool", m_cellPositionsHCalBarrelNoSegTool, - "Handle for tool to retrieve cell positions in HCal Barrel without DD4hep segmentation"); - declareProperty("positionsHCalExtBarrelTool", m_cellPositionsHCalExtBarrelTool, - "Handle for tool to retrieve cell positions in HCal ext Barrel"); - //declareProperty("positionsEMECTool", m_cellPositionsEMECTool, "Handle for tool to retrieve cell positions in EMEC"); - //declareProperty("positionsHECTool", m_cellPositionsHECTool, "Handle for tool to retrieve cell positions in HEC"); - //declareProperty("positionsEMFwdTool", m_cellPositionsEMFwdTool, "Handle for tool to retrieve cell positions EM Fwd"); - //declareProperty("positionsHFwdTool", m_cellPositionsHFwdTool, "Handle for tool to retrieve cell positions Had Fwd"); - declareProperty("clusters", m_clusterCollection, "Handle for calo clusters (output collection)"); - declareProperty("clusterCells", m_clusterCellsCollection, "Handle for clusters (output collection)"); + declareProperty("cells", m_cells, + "Handle for cells from ECal Barrel (input collection)"); + declareProperty("clusters", m_clusterCollection, + "Handle for calo clusters (output collection)"); + declareProperty("clusterCells", m_clusterCellsCollection, + "Handle for clusters (output collection)"); + declareProperty("noiseTool", m_noiseTool, + "Handle for the cells noise tool"); + declareProperty("neigboursTool", m_neighboursTool, + "Handle for tool to retrieve cell neighbors"); } StatusCode CaloTopoClusterFCCee::initialize() { - if (GaudiAlgorithm::initialize().isFailure()) return StatusCode::FAILURE; + if (GaudiAlgorithm::initialize().isFailure()) { + return StatusCode::FAILURE; + } + + // Check geometry service m_geoSvc = service("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_inputTool.retrieve()) { - error() << "Unable to retrieve the topo cluster input tool!!!" << endmsg; - return StatusCode::FAILURE; - } + + // Check the neighbours tool if (!m_neighboursTool.retrieve()) { error() << "Unable to retrieve the cells neighbours tool!!!" << endmsg; return StatusCode::FAILURE; } + + // Check the noise tool if (!m_noiseTool.retrieve()) { error() << "Unable to retrieve the cells noise tool!!!" << endmsg; return StatusCode::FAILURE; } - // Check if cell position ECal Barrel tool available - if (!m_cellPositionsECalBarrelTool.retrieve()) { - error() << "Unable to retrieve ECal Barrel cell positions tool!!!" << endmsg; - return StatusCode::FAILURE; - } - // Check if cell position HCal Barrel tool available - if (!m_cellPositionsHCalBarrelTool.retrieve()) { - error() << "Unable to retrieve HCal Barrel cell positions tool!!!" << endmsg; - if (!m_cellPositionsHCalBarrelNoSegTool.retrieve()) { - error() << "Also unable to retrieve HCal Barrel no segmentation cell positions tool!!!" << endmsg; - return StatusCode::FAILURE; - } - } - m_decoder_ecal = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder(); - m_index_layer_ecal = m_decoder_ecal->index("layer"); + // Setup system decoder + m_decoder = new dd4hep::DDSegmentation::BitFieldCoder(m_systemEncoding); + m_indexSystem = m_decoder->index("system"); + + // Setup ECal Barrel decoder + if (!m_readoutName.value().empty()) { + m_decoder_ecal = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder(); + m_index_layer_ecal = m_decoder_ecal->index("layer"); + } return StatusCode::SUCCESS; } StatusCode CaloTopoClusterFCCee::execute() { - - //std::unordered_map allCells; // transform it into a member variable to avoid - // recreating such a huge map at every event. Only update the energy values - std::vector> firstSeeds; - - // get input cell map from input tool - StatusCode sc_prepareCellMap = m_inputTool->cellIDMap(m_allCells); - if (sc_prepareCellMap.isFailure()) { - error() << "Unable to create cell map!" << endmsg; - return StatusCode::FAILURE; - } - debug() << "Active Cells : " << m_allCells.size() << endmsg; + // Get the input collection with calorimeter cells + const edm4hep::CalorimeterHitCollection* inCells = m_cells.get(); + debug() << "Input calo cell collection size: " << inCells->size() << endmsg; - // on first event, create cache - if (m_min_noise.empty()) { - createCache(m_allCells); + // On first event, create cache + if (m_min_noise.empty() && m_decoder_ecal) { + createCache(inCells); } - - // Create output collections - auto edmClusters = m_clusterCollection.createAndPut(); - std::unique_ptr edmClusterCells(new edm4hep::CalorimeterHitCollection()); - - // Finds seeds - CaloTopoClusterFCCee::findingSeeds(m_allCells, m_seedSigma, firstSeeds); - debug() << "Number of seeds found : " << firstSeeds.size() << endmsg; - // decending order of seeds - std::sort(firstSeeds.begin(), firstSeeds.end(), - [](const std::pair& lhs, const std::pair& rhs) { - return lhs.second < rhs.second; - }); - - std::map>> preClusterCollection; - StatusCode sc_buildProtoClusters = CaloTopoClusterFCCee::buildingProtoCluster(m_neighbourSigma, - m_lastNeighbourSigma, - firstSeeds, m_allCells, - preClusterCollection); - if (sc_buildProtoClusters.isFailure()) { - error() << "Unable to build the protoclusters!" << endmsg; - return StatusCode::FAILURE; + // Create output collections + edm4hep::ClusterCollection* outClusters = m_clusterCollection.createAndPut(); + edm4hep::CalorimeterHitCollection* outClusterCells = m_clusterCellsCollection.createAndPut(); + + // Find seeds + edm4hep::CalorimeterHitCollection seedCells = findSeeds(inCells); + debug() << "Number of seeds found : " << seedCells.size() << endmsg; + + // Build protoclusters (find neighbouring cells) + std::map protoClusters; + { + StatusCode sc = buildProtoClusters(seedCells, + inCells, + protoClusters); + if (sc.isFailure()) { + error() << "Unable to build the protoclusters!" << endmsg; + return StatusCode::FAILURE; + } } - // Build Clusters in edm - debug() << "Building " << preClusterCollection.size() << " cluster." << endmsg; + // Build clusters + debug() << "Building " << protoClusters.size() << " clusters." << endmsg; double checkTotEnergy = 0.; int clusterWithMixedCells = 0; - for (auto i : preClusterCollection) { + for (const auto& protoCluster: protoClusters) { edm4hep::MutableCluster cluster; //auto& clusterCore = cluster.core(); - double posX = 0.; - double posY = 0.; - double posZ = 0.; - double energy = 0.; + double clusterPosX = 0.; + double clusterPosY = 0.; + double clusterPosZ = 0.; + double clusterEnergy = 0.; double deltaR = 0.; - std::vector posPhi (i.second.size()); - std::vector posTheta (i.second.size()); - std::vector vecEnergy (i.second.size()); - double sumPhi = 0.; - double sumTheta = 0.; - std::map system; - - for (auto pair : i.second) { - dd4hep::DDSegmentation::CellID cID = pair.first; - // auto cellID = pair.first; - // get CalorimeterHit by cellID - auto newCell = edmClusterCells->create(); - newCell.setEnergy(m_allCells[cID]); - newCell.setCellID(cID); - newCell.setType(pair.second); - energy += newCell.getEnergy(); - - // get cell position by cellID + std::vector cellPosPhi(protoCluster.second.size(), 0); + std::vector cellPosTheta(protoCluster.second.size(), 0); + std::vector cellEnergy(protoCluster.second.size(), 0); + double sumCellPhi = 0.; + double sumCellTheta = 0.; + std::map system; + + debug() << "Building cluster with ID: " << protoCluster.first << endmsg; + + for (const auto& protoCell: protoCluster.second) { + auto cell = protoCell.clone(); + clusterEnergy += cell.getEnergy(); // identify calo system - auto systemId = m_decoder->get(cID, m_index_system); + auto systemId = m_decoder->get(cell.getCellID(), m_indexSystem); system[int(systemId)]++; - dd4hep::Position posCell; - if (systemId == 4) // ECAL BARREL system id - posCell = m_cellPositionsECalBarrelTool->xyzPosition(cID); - else if (systemId == 8){ // HCAL BARREL system id - if (m_noSegmentationHCalUsed) - posCell = m_cellPositionsHCalBarrelNoSegTool->xyzPosition(cID); - else - posCell = m_cellPositionsHCalBarrelTool->xyzPosition(cID); - } - //else if (systemId == 9) // HCAL EXT BARREL system id - // posCell = m_cellPositionsHCalExtBarrelTool->xyzPosition(cID); - //else if (systemId == 6) // EMEC system id - // posCell = m_cellPositionsEMECTool->xyzPosition(cID); - //else if (systemId == 7) // HEC system id - // posCell = m_cellPositionsHECTool->xyzPosition(cID); - //else if (systemId == 10) // EMFWD system id - // posCell = m_cellPositionsEMFwdTool->xyzPosition(cID); - //else if (systemId == 11) // HFWD system id - // posCell = m_cellPositionsHFwdTool->xyzPosition(cID); - else - warning() << "No cell positions tool found for system id " << systemId << ". " << endmsg; - - posX += posCell.X() * newCell.getEnergy(); - posY += posCell.Y() * newCell.getEnergy(); - posZ += posCell.Z() * newCell.getEnergy(); - posPhi.push_back(posCell.Phi()); - posTheta.push_back(posCell.Theta()); - vecEnergy.push_back(newCell.getEnergy()); - sumPhi += posCell.Phi() * newCell.getEnergy(); - sumTheta += posCell.Theta() * newCell.getEnergy(); - - cluster.addToHits(newCell); - auto er = m_allCells.erase(cID); - - if (er!=1) - info() << "Problem in erasing cell ID from map." << endmsg; + auto cellPos = dd4hep::Position(cell.getPosition().x, + cell.getPosition().y, + cell.getPosition().z); + + clusterPosX += cell.getPosition().x * cell.getEnergy(); + clusterPosY += cell.getPosition().y * cell.getEnergy(); + clusterPosZ += cell.getPosition().z * cell.getEnergy(); + cellPosPhi.push_back(cellPos.Phi()); + cellPosTheta.push_back(cellPos.Theta()); + cellEnergy.push_back(cell.getEnergy()); + sumCellPhi += cellPos.Phi() * cell.getEnergy(); + sumCellTheta += cellPos.Theta() * cell.getEnergy(); + + cluster.addToHits(cell); + outClusterCells->push_back(cell); } - cluster.setEnergy(energy); - cluster.setPosition(edm4hep::Vector3f(posX / energy, posY / energy, posZ / energy)); + cluster.setEnergy(clusterEnergy); + cluster.setPosition(edm4hep::Vector3f(clusterPosX / clusterEnergy, + clusterPosY / clusterEnergy, + clusterPosZ / clusterEnergy)); // store deltaR of cluster in time for the moment.. - sumPhi = sumPhi / energy; - sumTheta = sumTheta / energy; - int counter = 0; - for (auto entryTheta : posTheta){ - deltaR += sqrt(pow(entryTheta-sumTheta,2) + pow(posPhi[counter]-sumPhi,2)) * vecEnergy[counter]; - counter++; + sumCellPhi = sumCellPhi / clusterEnergy; + sumCellTheta = sumCellTheta / clusterEnergy; + for (size_t i = 0; i < cellEnergy.size(); ++i) { + deltaR += std::sqrt(std::pow(cellPosTheta[i] - sumCellTheta, 2) + + std::pow(cellPosPhi[i] - sumCellPhi, 2)) + * cellEnergy[i]; } - cluster.addToShapeParameters(deltaR / energy); + cluster.addToShapeParameters(deltaR / clusterEnergy); verbose() << "Cluster energy: " << cluster.getEnergy() << endmsg; checkTotEnergy += cluster.getEnergy(); - edmClusters->push_back(cluster); + outClusters->push_back(cluster); if (system.size() > 1) clusterWithMixedCells++; - posPhi.clear(); - posTheta.clear(); - vecEnergy.clear(); - + cellPosPhi.clear(); + cellPosTheta.clear(); + cellEnergy.clear(); } - m_clusterCellsCollection.put(std::move(edmClusterCells)); - debug() << "Number of clusters with cells in E and HCal: " << clusterWithMixedCells << endmsg; - debug() << "Total energy of clusters: " << checkTotEnergy << endmsg; - debug() << "Leftover cells : " << m_allCells.size() << endmsg; + debug() << "Number of clusters with cells in E and HCal: " + << clusterWithMixedCells << endmsg; + debug() << "Total energy of clusters: " + << checkTotEnergy << endmsg; + debug() << "Leftover cells : " + << inCells->size() - outClusterCells->size() << endmsg; + return StatusCode::SUCCESS; } -void CaloTopoClusterFCCee::findingSeeds(const std::unordered_map& aCells, - int aNumSigma, - std::vector>& aSeeds) { - for (const auto& cell : aCells) { + +StatusCode CaloTopoClusterFCCee::finalize() { return GaudiAlgorithm::finalize(); } + + +edm4hep::CalorimeterHitCollection CaloTopoClusterFCCee::findSeeds( + const edm4hep::CalorimeterHitCollection* allCells) { + + std::vector seedCellsVec; + + for (const auto& cell : *allCells) { // retrieve the noise const and offset assigned to cell // first try to use the cache - int system = m_decoder->get(cell.first, m_index_system); - if (system == 4) { //ECal barrel - int layer = m_decoder_ecal->get(cell.first, m_index_layer_ecal); + int system = m_decoder->get(cell.getCellID(), m_indexSystem); + if (system == 4 && m_decoder_ecal) { // ECal barrel + int layer = m_decoder_ecal->get(cell.getCellID(), m_index_layer_ecal); - double min_threshold = m_min_offset[layer] + m_min_noise[layer] * aNumSigma; + double min_threshold = m_min_offset[layer] + + m_min_noise[layer] * m_seedSigma; debug() << "m_min_offset[layer] = " << m_min_offset[layer] << endmsg; debug() << "m_min_noise[layer] = " << m_min_noise[layer] << endmsg; - debug() << "aNumSigma = " << aNumSigma << endmsg; + debug() << "aNumSigma = " << m_seedSigma << endmsg; debug() << "min_threshold = " << min_threshold << endmsg; - debug() << "abs(cell.second) = " << abs(cell.second) << endmsg; + debug() << "abs(cell.second) = " << std::fabs(cell.getEnergy()) << endmsg; - if (abs(cell.second) < min_threshold) { - // if we are below the minimum threshold for the full layer, no need to retrieve the exact value + if (std::fabs(cell.getEnergy()) < min_threshold) { + // if we are below the minimum threshold for the full layer, no need to + // retrieve the exact value continue; } } - // we are above the minimum threshold of the layer. Let's see if we are above the threshold for this cell. - double threshold = m_noiseTool->noiseOffset(cell.first) + ( m_noiseTool->noiseRMS(cell.first) * aNumSigma); - if (msgLevel() <= MSG::VERBOSE){ - debug() << "noise offset = " << m_noiseTool->noiseOffset(cell.first) << "GeV " << endmsg; - debug() << "noise rms = " << m_noiseTool->noiseRMS(cell.first) << "GeV " << endmsg; + // we are above the minimum threshold of the layer. Let's see if we are + // above the threshold for this cell. + double threshold = m_noiseTool->noiseOffset(cell.getCellID()) + + (m_noiseTool->noiseRMS(cell.getCellID()) * m_seedSigma); + if (msgLevel() <= MSG::VERBOSE) { + debug() << "noise offset = " << m_noiseTool->noiseOffset(cell.getCellID()) << "GeV " << endmsg; + debug() << "noise rms = " << m_noiseTool->noiseRMS(cell.getCellID()) << "GeV " << endmsg; debug() << "seed threshold = " << threshold << "GeV " << endmsg; debug() << "======================================" << endmsg; } - if (abs(cell.second) > threshold) { - aSeeds.emplace_back(cell.first, cell.second); + if (std::fabs(cell.getEnergy()) > threshold) { + seedCellsVec.emplace_back(cell); } } + + // descending order of seeds + std::sort(seedCellsVec.begin(), seedCellsVec.end(), + [](const auto& lhs, const auto& rhs) { + return lhs.getEnergy() < rhs.getEnergy(); + }); + + edm4hep::CalorimeterHitCollection seedCells; + seedCells.setSubsetCollection(); + for (const auto& cell: seedCellsVec) { + seedCells.push_back(cell); + } + + return seedCells; } -StatusCode CaloTopoClusterFCCee::buildingProtoCluster( - int aNumSigma, - int aLastNumSigma, - std::vector>& aSeeds, - const std::unordered_map& aCells, - std::map>>& aPreClusterCollection) { - // Map of cellIDs to clusterIds - std::map clusterOfCell; + +StatusCode CaloTopoClusterFCCee::buildProtoClusters( + const edm4hep::CalorimeterHitCollection& seedCells, + const edm4hep::CalorimeterHitCollection* allCells, + std::map& protoClusters) { + + verbose() << "Initial number of seeds to loop over: " << seedCells.size() + << endmsg; + + std::map allCellsMap; + for (const auto& cell: *allCells) { + allCellsMap.emplace(cell.getCellID(), cell); + } + std::map alreadyUsedCells; // Loop over every seed in Calo to create first cluster - uint iSeeds = 0; - verbose() << "seeds to loop over : " << aSeeds.size() << endmsg; - for (auto& itSeed : aSeeds) { - iSeeds++; - verbose() << "Seed num: " << iSeeds << endmsg; - auto seedId = itSeed.first; - auto cellInCluster = clusterOfCell.find(seedId); - if (cellInCluster != clusterOfCell.end()) { + uint32_t seedCounter = 0; + for (const auto& seedCell: seedCells) { + seedCounter++; + verbose() << "Looking at seed: " << seedCounter << endmsg; + auto seedId = seedCell.getCellID(); + auto cellInCluster = alreadyUsedCells.find(seedId); + if (cellInCluster != alreadyUsedCells.end()) { verbose() << "Seed is already assigned to another cluster!" << endmsg; continue; - } else { - // new cluster starts with seed - // set cell Bits to 1 for seed cell - aPreClusterCollection[iSeeds].push_back(std::make_pair(seedId, 1)); - uint clusterId = iSeeds; - clusterOfCell[seedId] = clusterId; - - std::vector>> vecNextNeighbours(1); - vecNextNeighbours.reserve(100); - vecNextNeighbours[0] = CaloTopoClusterFCCee::searchForNeighbours(seedId, clusterId, aNumSigma, aCells, clusterOfCell, - aPreClusterCollection, true); - // first loop over seeds neighbours - verbose() << "Found " << vecNextNeighbours[0].size() << " neighbours.." << endmsg; - int it = 0; - while (vecNextNeighbours[it].size() > 0) { - it++; - vecNextNeighbours.emplace_back(std::vector>{}); - for (auto& id : vecNextNeighbours[it - 1]) { - if (id.first == 0){ - error() << "Building of cluster is stopped due to missing id in neighbours map." << endmsg; - return StatusCode::FAILURE; - } - verbose() << "Next neighbours assigned to clusterId : " << clusterId << endmsg; - auto vec = CaloTopoClusterFCCee::searchForNeighbours(id.first, clusterId, aNumSigma, aCells, clusterOfCell, - aPreClusterCollection, true); - vecNextNeighbours[it].insert(vecNextNeighbours[it].end(), vec.begin(), vec.end()); + } + + uint32_t clusterId = seedCounter; + // new cluster starts with seed + // set cell type to 1 for seed cell + edm4hep::MutableCalorimeterHit clusteredCell = seedCell.clone(); + clusteredCell.setType(1); + protoClusters[clusterId].push_back(clusteredCell); + alreadyUsedCells[seedId] = clusterId; + + std::vector>> nextNeighbours(100); + nextNeighbours[0] = searchForNeighbours(seedId, + clusterId, + m_neighbourSigma, + allCellsMap, + alreadyUsedCells, + protoClusters, + true); + + // first loop over seeds neighbours + verbose() << "Found " << nextNeighbours[0].size() << " neighbours.." + << endmsg; + int it = 0; + while (nextNeighbours[it].size() > 0) { + it++; + nextNeighbours.emplace_back(std::vector>{}); + for (auto& id : nextNeighbours[it - 1]) { + if (id.first == 0) { + error() << "Building of cluster is stopped due to missing cell ID " + "in neighbours map!" << endmsg; + return StatusCode::FAILURE; } - verbose() << "Found " << vecNextNeighbours[it].size() << " more neighbours.." << endmsg; + verbose() << "Next neighbours assigned to cluster ID: " << clusterId + << endmsg; + auto additionalNeighbours = searchForNeighbours(id.first, + clusterId, + m_neighbourSigma, + allCellsMap, + alreadyUsedCells, + protoClusters, + true); + nextNeighbours[it].insert(nextNeighbours[it].end(), + additionalNeighbours.begin(), + additionalNeighbours.end()); } - // last try with different condition on neighbours - if (vecNextNeighbours[it].size() == 0) { - auto clusteredCells = aPreClusterCollection[clusterId]; - // loop over all clustered cells - for (auto& id : clusteredCells) { - if (id.second <= 2){ - verbose() << "Add neighbours of " << id.first << " in last round with thr = " << aLastNumSigma << " x sigma." << endmsg; - auto lastNeighours = CaloTopoClusterFCCee::searchForNeighbours(id.first, clusterId, aLastNumSigma, aCells, clusterOfCell, - aPreClusterCollection, false); - } - } + verbose() << "Found " << nextNeighbours[it].size() + << " more neighbours.." << endmsg; + } + + // last try with different condition on neighbours + if (nextNeighbours[it].size() == 0) { + // loop over all clustered cells + for (const auto& cell: protoClusters[clusterId]) { + if (cell.getType() <= 2) { + verbose() << "Add neighbours of " << cell.getCellID() + << " in last round with thr = " + << m_lastNeighbourSigma.value() + << " x sigma." << endmsg; + auto lastNeighours = searchForNeighbours(cell.getCellID(), + clusterId, + m_lastNeighbourSigma, + allCellsMap, + alreadyUsedCells, + protoClusters, + false); + } } } } + return StatusCode::SUCCESS; } -std::vector > -CaloTopoClusterFCCee::searchForNeighbours(const uint64_t aCellId, - uint& aClusterID, - int aNumSigma, - const std::unordered_map& aCells, - std::map& aClusterOfCell, - std::map>>& aPreClusterCollection, - bool aAllowClusterMerge) { - // Fill vector to be returned, next cell ids and cluster id for which neighbours are found - std::vector> addedNeighbourIds; +std::vector> +CaloTopoClusterFCCee::searchForNeighbours ( + const uint64_t aCellId, + uint32_t& aClusterID, + const int aNumSigma, + std::map& allCellsMap, + std::map& alreadyUsedCells, + std::map& protoClusters, + const bool aAllowClusterMerge) { + // Fill vector to be returned, next cell ids and cluster id for which + // neighbours are found + std::vector> additionalNeighbours; + // Retrieve cellIDs of neighbours auto neighboursVec = m_neighboursTool->neighbours(aCellId); if (neighboursVec.size() == 0) { error() << "No neighbours for cellID found! " << endmsg; error() << "to cellID : " << aCellId << endmsg; - error() << "in system: " << m_decoder->get(aCellId, m_index_system) << endmsg; - addedNeighbourIds.resize(0); - addedNeighbourIds.push_back(std::make_pair(0, 0)); - } else { - - verbose() << "For cluster: " << aClusterID << endmsg; - // loop over neighbours - for (auto& itr : neighboursVec) { - auto neighbourID = itr; - // Find the neighbour in the Calo cells list - auto itAllCells = aCells.find(neighbourID); - auto itAllUsedCells = aClusterOfCell.find(neighbourID); - - // If cell is hit.. and is not assigned to a cluster - if (itAllCells != aCells.end() && itAllUsedCells == aClusterOfCell.end()) { - verbose() << "Found neighbour with CellID: " << neighbourID << endmsg; - auto neighbouringCellEnergy = itAllCells->second; - bool addNeighbour = false; - int cellType = 2; - // retrieve the cell noise level [GeV] - // - // first try to use the cache - bool is_below = false; - int system = m_decoder->get(neighbourID, m_index_system); - if (system == 4) { //ECal barrel - int layer = m_decoder_ecal->get(neighbourID, m_index_layer_ecal); - double min_threshold = m_min_offset[layer] + m_min_noise[layer] * aNumSigma; - if (abs(neighbouringCellEnergy) < min_threshold) { - // if we are below the minimum threshold for the full layer, no need to retrieve the exact value - is_below = true; - } + error() << "in system: " << m_decoder->get(aCellId, m_indexSystem) << endmsg; + additionalNeighbours.resize(0); + additionalNeighbours.push_back(std::make_pair(0, 0)); + return additionalNeighbours; } - if (is_below) { - addNeighbour = false; - } - else { - double thr = m_noiseTool->noiseOffset(neighbourID) + (aNumSigma * m_noiseTool->noiseRMS(neighbourID)); - if (abs(neighbouringCellEnergy) > thr) - addNeighbour = true; - else - addNeighbour = false; - } - // give cell type according to threshold - if (aNumSigma == m_lastNeighbourSigma){ - cellType = 3; - } - // if threshold is 0, collect the cell independent on its energy - if (aNumSigma == 0){ - addNeighbour = true; - } - // if neighbour is validated - if (addNeighbour) { - // retrieve the cell - // add neighbour to cells for cluster - aPreClusterCollection[aClusterID].push_back(std::make_pair(neighbourID, cellType)); - aClusterOfCell[neighbourID] = aClusterID; - addedNeighbourIds.push_back(std::make_pair(neighbourID, aClusterID)); + verbose() << "For cluster: " << aClusterID << endmsg; + // loop over neighbours + for (const auto& neighbourID : neighboursVec) { + // Find the neighbour in the Calo cells list + auto itAllCells = allCellsMap.find(neighbourID); + auto itAllUsedCells = alreadyUsedCells.find(neighbourID); + + // If cell is hit.. and is not assigned to a cluster + if (itAllCells != allCellsMap.end() && itAllUsedCells == alreadyUsedCells.end()) { + verbose() << "Found neighbour with CellID: " << neighbourID << endmsg; + auto neighbouringCellEnergy = allCellsMap[neighbourID].getEnergy(); + bool addNeighbour = false; + int cellType = 2; + // retrieve the cell noise level [GeV] + // + // first try to use the cache + bool isBelow = false; + int system = m_decoder->get(neighbourID, m_indexSystem); + if (system == 4 && m_decoder_ecal) { // ECal barrel + int layer = m_decoder_ecal->get(neighbourID, m_index_layer_ecal); + double min_threshold = m_min_offset[layer] + + m_min_noise[layer] * aNumSigma; + if (std::fabs(neighbouringCellEnergy) < min_threshold) { + // if we are below the minimum threshold for the full layer, no + // need to retrieve the exact value + isBelow = true; } } - // If cell is hit.. but is assigned to another cluster - else if (itAllUsedCells != aClusterOfCell.end() && itAllUsedCells->second != aClusterID && aAllowClusterMerge) { - uint clusterIDToMerge = itAllUsedCells->second; - if (msgLevel() <= MSG::VERBOSE){ - verbose() << "This neighbour was found in cluster " << clusterIDToMerge << ", cluster " << aClusterID - << " will be merged!" << endmsg; - verbose() << "Assigning all cells ( " << aPreClusterCollection.find(aClusterID)->second.size() << " ) to Cluster " - << clusterIDToMerge << " with ( " << aPreClusterCollection.find(clusterIDToMerge)->second.size() - << " ). " << endmsg; - } - // Fill all cells into cluster, and assigned cells to new cluster - aClusterOfCell[neighbourID] = clusterIDToMerge; - for (auto& i : aPreClusterCollection.find(aClusterID)->second) { - aClusterOfCell[i.first] = clusterIDToMerge; - bool found = false; - // make sure that already assigned cells are not added - for (auto& j : aPreClusterCollection[clusterIDToMerge]) { - if (j.first == i.first) found = true; - } - if (!found) { - aPreClusterCollection[clusterIDToMerge].push_back(std::make_pair(i.first, i.second)); - } + + if (isBelow) { + addNeighbour = false; + } + else { + double thr = m_noiseTool->noiseOffset(neighbourID) + + (aNumSigma * m_noiseTool->noiseRMS(neighbourID)); + if (std::fabs(neighbouringCellEnergy) > thr) + addNeighbour = true; + else + addNeighbour = false; + } + // give cell type according to threshold + if (aNumSigma == m_lastNeighbourSigma){ + cellType = 3; + } + // if threshold is 0, collect the cell independent on its energy + if (aNumSigma == 0){ + addNeighbour = true; + } + // if neighbour is validated + if (addNeighbour) { + // retrieve the cell + // add neighbour to cells for cluster + edm4hep::MutableCalorimeterHit clusteredCell = + allCellsMap[neighbourID].clone(); + clusteredCell.setType(cellType); + protoClusters[aClusterID].push_back(clusteredCell); + alreadyUsedCells[neighbourID] = aClusterID; + additionalNeighbours.push_back(std::make_pair(neighbourID, aClusterID)); + } + } + // If cell is hit.. but is assigned to another cluster + else if (itAllUsedCells != alreadyUsedCells.end() && itAllUsedCells->second != aClusterID && aAllowClusterMerge) { + uint32_t clusterIDToMergeTo = itAllUsedCells->second; + if (msgLevel() <= MSG::VERBOSE){ + verbose() << "This neighbour was found in cluster " << clusterIDToMergeTo + << ", cluster " << aClusterID + << " will be merged!" << endmsg; + verbose() << "Assigning all cells ( " + << protoClusters[aClusterID].size() << " ) to Cluster " + << clusterIDToMergeTo << " with ( " + << protoClusters[clusterIDToMergeTo].size() + << " ). " << endmsg; + } + // Fill all cells into cluster, and assigned cells to new cluster + alreadyUsedCells[neighbourID] = clusterIDToMergeTo; + for (const auto& cell : protoClusters[aClusterID]) { + alreadyUsedCells[cell.getCellID()] = clusterIDToMergeTo; + // make sure that already assigned cells are not added + if (cellIdInColl(cell.getCellID(), protoClusters[clusterIDToMergeTo])) { + continue; } - aPreClusterCollection.erase(aClusterID); - // changed clusterId -> if more neighbours are found, correct assignment - verbose() << "Cluster Id changed to " << clusterIDToMerge << endmsg; - aClusterID = clusterIDToMerge; - // found neighbour for next search - addedNeighbourIds.push_back(std::make_pair(neighbourID, aClusterID)); - // end loop to ensure correct cluster assignment - break; + protoClusters[clusterIDToMergeTo].push_back(cell.clone()); } + protoClusters.erase(aClusterID); + // changed clusterId -> if more neighbours are found, correct assignment + verbose() << "Cluster Id changed to " << clusterIDToMergeTo << endmsg; + aClusterID = clusterIDToMergeTo; + // found neighbour for next search + additionalNeighbours.push_back(std::make_pair(neighbourID, aClusterID)); + // end loop to ensure correct cluster assignment + break; } } - return addedNeighbourIds; + + return additionalNeighbours; } -StatusCode CaloTopoClusterFCCee::finalize() { return GaudiAlgorithm::finalize(); } -void CaloTopoClusterFCCee::createCache(const std::unordered_map& aCells) { +void CaloTopoClusterFCCee::createCache( + const edm4hep::CalorimeterHitCollection* aCells) { // cache the minimum offset and noise per layer for faster lookups down the chain. std::unordered_map> offsets; @@ -454,17 +483,17 @@ void CaloTopoClusterFCCee::createCache(const std::unordered_map layers; // fill all noises and offsets values - for (const auto& cell : aCells) { - int system = m_decoder->get(cell.first, m_index_system); - if (system == 4) { //ECal barrel - int layer = m_decoder_ecal->get(cell.first, m_index_layer_ecal); + for (const auto& cell : *aCells) { + int system = m_decoder->get(cell.getCellID(), m_indexSystem); + if (system == 4 && m_decoder_ecal) { //ECal barrel + int layer = m_decoder_ecal->get(cell.getCellID(), m_index_layer_ecal); if (layers.find(layer) == layers.end()) { offsets[layer] = std::vector{}; noises[layer] = std::vector{}; layers.insert(layer); } - offsets[layer].push_back(m_noiseTool->noiseOffset(cell.first)); - noises[layer].push_back(m_noiseTool->noiseRMS(cell.first)); + offsets[layer].push_back(m_noiseTool->noiseOffset(cell.getCellID())); + noises[layer].push_back(m_noiseTool->noiseRMS(cell.getCellID())); } } @@ -479,6 +508,17 @@ void CaloTopoClusterFCCee::createCache(const std::unordered_map +#include +#include +#include +#include #include class IGeoSvc; @@ -50,28 +53,25 @@ class CaloTopoClusterFCCee : public GaudiAlgorithm { StatusCode initialize(); - /** Find cells with a signal to noise ratio > nSigma. - * @param[in] aCells, parse the map of all cells. - * @param[in] aNumSigma, the signal to noise ratio that the cell has to exceed to become seed. - * @param[in] aSeeds, the vector of seed cell ids anf their energy to build proto-clusters. + /** Find cells with a signal to noise ratio > seedSigma. + * @param[in] allCells, collection of all cells. + * return collection of seed cells. */ - virtual void findingSeeds(const std::unordered_map& aCells, int aNumSigma, - std::vector>& aSeeds); + edm4hep::CalorimeterHitCollection findSeeds( + const edm4hep::CalorimeterHitCollection* allCells); /** Building proto-clusters from the found seeds. * First the function initialises a cluster in the preClusterCollection for the seed cells, * then it calls the CaloTopoClusterFCCee::searchForNeighbours function to retrieve the vector of next cellIDs to add and loop over to find neighbours. * The iteration of search for neighbours is continued until no more neihgbours are found. Then a last round of adding neighbouring cells to the cluster is run where the parameter lastNeighbourSigma is applied. - * @param[in] aNumSigma, signal to noise ratio the neighbouring cell has to pass to be added to cluster. * @param[in] aSeeds, vector of seeding cells. * @param[in] aCells, map of all cells. * @param[in] aPreClusterCollection, map that is filled with clusterID pointing to the associated cells, in a pair of cellID and cellType. */ - StatusCode buildingProtoCluster(int aNumSigma, - int aLastNumSigma, - std::vector>& aSeeds, - const std::unordered_map& aCells, - std::map>>& aPreClusterCollection); + StatusCode buildProtoClusters( + const edm4hep::CalorimeterHitCollection& seedCells, + const edm4hep::CalorimeterHitCollection* allCells, + std::map& protoClusters); /** Search for neighbours and add them to preClusterCollection * The @@ -84,47 +84,40 @@ class CaloTopoClusterFCCee : public GaudiAlgorithm { * @param[in] aAllowClusterMerge, bool to allow for clusters to be merged, set to false in case of last iteration in CaloTopoClusterFCCee::buildingProtoCluster. * return vector of pairs with cellID and energy of found neighbours. */ - std::vector> - searchForNeighbours(const uint64_t aCellId, uint& aClusterID, int aNumSigma, const std::unordered_map& aCells, - std::map& aClusterOfCell, - std::map>>& aPreClusterCollection, - bool aAllowClusterMerge); + std::vector> + searchForNeighbours( + const uint64_t aCellId, + uint32_t& aClusterID, + const int aNumSigma, + std::map& aCellsMap, + std::map& aClusterOfCell, + std::map& protoClusters, + const bool aAllowClusterMerge); StatusCode execute(); StatusCode finalize(); private: - // Cluster collection - DataHandle m_clusterCollection{"calo/clusters", Gaudi::DataHandle::Writer, this}; - // Cluster cells in collection - DataHandle m_clusterCellsCollection{"calo/clusterCells", Gaudi::DataHandle::Writer, this}; /// Pointer to the geometry service SmartIF m_geoSvc; - /// Handle for the input tool - ToolHandle m_inputTool{"TopoClusterInput", this}; + + /// Handle for electromagnetic barrel cells (input collection) + DataHandle m_cells{ + "cells", Gaudi::DataHandle::Reader, this}; + // Cluster collection (output) + DataHandle m_clusterCollection{ + "clusters", Gaudi::DataHandle::Writer, this}; + // Cluster cells in collection (output) + DataHandle m_clusterCellsCollection{ + "clusterCells", Gaudi::DataHandle::Writer, this}; + /// Handle for the cells noise tool ToolHandle m_noiseTool{"TopoCaloNoisyCells", this}; - /// Handle for neighbours tool + /// Handle for neighbors tool ToolHandle m_neighboursTool{"TopoCaloNeighbours", this}; - /// Handle for tool to get positions in ECal Barrel - ToolHandle m_cellPositionsECalBarrelTool{"CellPositionsECalBarrelTool", this}; - /// Handle for tool to get positions in HCal Barrel - ToolHandle m_cellPositionsHCalBarrelNoSegTool{"CellPositionsHCalBarrelNoSegTool", this}; - /// Handle for tool to get positions in HCal Barrel - ToolHandle m_cellPositionsHCalBarrelTool{"CellPositionsHCalBarrelTool", this}; - /// Handle for tool to get positions in HCal Barrel and Ext Barrel, no Segmentation - ToolHandle m_cellPositionsHCalExtBarrelTool{"CellPositionsHCalBarrelNoSegTool", this}; - ///// Handle for tool to get positions in Calo Discs - //ToolHandle m_cellPositionsEMECTool{"CellPositionsCaloDiscsTool", this}; - ///// Handle for tool to get positions in Calo Discs - //ToolHandle m_cellPositionsHECTool{"CellPositionsCaloDiscsTool", this}; - ///// Handle for tool to get positions in Calo Discs - //ToolHandle m_cellPositionsEMFwdTool{"CellPositionsCaloDiscsTool", this}; - ///// Handle for tool to get positions in Calo Discs - //ToolHandle m_cellPositionsHFwdTool{"CellPositionsCaloDiscsTool", this}; - - /// no segmentation used in HCal + + /// No segmentation used in HCal Gaudi::Property m_noSegmentationHCalUsed{this, "noSegmentationHCal", true, "HCal Barrel readout without DD4hep theta-module segmentation used."}; /// Seed threshold in sigma Gaudi::Property m_seedSigma{this, "seedSigma", 4, "number of sigma in noise threshold"}; @@ -133,24 +126,30 @@ class CaloTopoClusterFCCee : public GaudiAlgorithm { /// Last neighbour threshold in sigma Gaudi::Property m_lastNeighbourSigma{this, "lastNeighbourSigma", 0, "number of sigma in noise threshold"}; /// Name of the electromagnetic calorimeter readout - Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelModuleThetaMerged"}; + // Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelModuleThetaMerged"}; + Gaudi::Property m_readoutName{this, "readoutName", ""}; - /// General decoder to encode the calorimeter sub-system to determine which positions tool to use - dd4hep::DDSegmentation::BitFieldCoder* m_decoder = new dd4hep::DDSegmentation::BitFieldCoder("system:4"); - int m_index_system = m_decoder->index("system"); + /// System encoding string + Gaudi::Property m_systemEncoding{ + this, "systemEncoding", "system:4", "System encoding string"}; + /// General decoder to encode the calorimeter sub-system to determine which + /// positions tool to use + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + int m_indexSystem; /// Decoder for ECal layers dd4hep::DDSegmentation::BitFieldCoder* m_decoder_ecal; int m_index_layer_ecal; - std::unordered_map m_allCells; - // minimum noise and offset per barrel ECal layer // this serves as a very small cache for fast lookups and avoid looking into the huge map for most of the cells. std::vector m_min_offset; std::vector m_min_noise; - void createCache(const std::unordered_map& aCells); + // Utility functions + void createCache(const edm4hep::CalorimeterHitCollection* aCells); + inline bool cellIdInColl(const uint64_t cellId, + const edm4hep::CalorimeterHitCollection& coll); }; #endif /* RECFCCEECALORIMETER_CALOTOPOCLUSTERFCCEE_H */