Skip to content

Commit

Permalink
Adjustments
Browse files Browse the repository at this point in the history
  • Loading branch information
@kjvbrt
kjvbrt committed Nov 9, 2023
1 parent c04343f commit 71106d5
Show file tree
Hide file tree
Showing 2 changed files with 60 additions and 33 deletions.
72 changes: 47 additions & 25 deletions RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,12 @@ StatusCode CaloTopoClusterFCCee::initialize() {
error() << "Unable to retrieve the cells noise tool!!!" << endmsg;
return StatusCode::FAILURE;
}

// Setup system decoder
m_decoder = new dd4hep::DDSegmentation::BitFieldCoder(m_systemEncoding);
m_indexSystem = m_decoder->index("system");

// Setup ECal Barrel decoder
m_decoder_ecal = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder();
m_index_layer_ecal = m_decoder_ecal->index("layer");

Expand Down Expand Up @@ -117,10 +123,11 @@ StatusCode CaloTopoClusterFCCee::execute() {

debug() << "Building cluster with ID: " << protoCluster.first << endmsg;

for (const auto& cell: protoCluster.second) {
for (const auto& protoCell: protoCluster.second) {
auto cell = protoCell.clone();
clusterEnergy += cell.getEnergy();
// identify calo system
auto systemId = m_decoder->get(cell.getCellID(), m_index_system);
auto systemId = m_decoder->get(cell.getCellID(), m_indexSystem);
system[int(systemId)]++;
auto cellPos = dd4hep::Position(cell.getPosition().x,
cell.getPosition().y,
Expand Down Expand Up @@ -168,7 +175,7 @@ StatusCode CaloTopoClusterFCCee::execute() {
debug() << "Total energy of clusters: "
<< checkTotEnergy << endmsg;
debug() << "Leftover cells : "
<< outClusterCells->size() - inCells->size() << endmsg;
<< inCells->size() - outClusterCells->size() << endmsg;

return StatusCode::SUCCESS;
}
Expand All @@ -185,7 +192,7 @@ edm4hep::CalorimeterHitCollection CaloTopoClusterFCCee::findSeeds(
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.getCellID(), m_index_system);
int system = m_decoder->get(cell.getCellID(), m_indexSystem);
if (system == 4) { // ECal barrel
int layer = m_decoder_ecal->get(cell.getCellID(), m_index_layer_ecal);

Expand All @@ -196,7 +203,7 @@ edm4hep::CalorimeterHitCollection CaloTopoClusterFCCee::findSeeds(
debug() << "m_min_noise[layer] = " << m_min_noise[layer] << endmsg;
debug() << "aNumSigma = " << m_seedSigma << endmsg;
debug() << "min_threshold = " << min_threshold << endmsg;
debug() << "abs(cell.second) = " << abs(cell.getEnergy()) << endmsg;
debug() << "abs(cell.second) = " << std::fabs(cell.getEnergy()) << endmsg;

if (std::fabs(cell.getEnergy()) < min_threshold) {
// if we are below the minimum threshold for the full layer, no need to
Expand Down Expand Up @@ -251,17 +258,18 @@ StatusCode CaloTopoClusterFCCee::buildProtoClusters(
std::map<uint64_t, uint32_t> alreadyUsedCells;

// Loop over every seed in Calo to create first cluster
uint32_t clusterId = 0;
uint32_t seedCounter = 0;
for (const auto& seedCell: seedCells) {
verbose() << "Clusters so far: " << clusterId << endmsg;
auto seedId = seedCell.getObjectID().index;
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;
}

clusterId++;
uint32_t clusterId = seedCounter;
// new cluster starts with seed
// set cell type to 1 for seed cell
edm4hep::MutableCalorimeterHit clusteredCell = seedCell.clone();
Expand Down Expand Up @@ -314,7 +322,8 @@ StatusCode CaloTopoClusterFCCee::buildProtoClusters(
for (const auto& cell: protoClusters[clusterId]) {
if (cell.getType() <= 2) {
verbose() << "Add neighbours of " << cell.getCellID()
<< " in last round with thr = " << m_lastNeighbourSigma
<< " in last round with thr = "
<< m_lastNeighbourSigma.value()
<< " x sigma." << endmsg;
auto lastNeighours = searchForNeighbours(cell.getCellID(),
clusterId,
Expand Down Expand Up @@ -349,7 +358,7 @@ CaloTopoClusterFCCee::searchForNeighbours (
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;
error() << "in system: " << m_decoder->get(aCellId, m_indexSystem) << endmsg;
additionalNeighbours.resize(0);
additionalNeighbours.push_back(std::make_pair(0, 0));
return additionalNeighbours;
Expand All @@ -372,7 +381,7 @@ CaloTopoClusterFCCee::searchForNeighbours (
//
// first try to use the cache
bool isBelow = false;
int system = m_decoder->get(neighbourID, m_index_system);
int system = m_decoder->get(neighbourID, m_indexSystem);
if (system == 4) { // ECal barrel
int layer = m_decoder_ecal->get(neighbourID, m_index_layer_ecal);
double min_threshold = m_min_offset[layer] +
Expand All @@ -390,7 +399,7 @@ CaloTopoClusterFCCee::searchForNeighbours (
else {
double thr = m_noiseTool->noiseOffset(neighbourID) +
(aNumSigma * m_noiseTool->noiseRMS(neighbourID));
if (abs(neighbouringCellEnergy) > thr)
if (std::fabs(neighbouringCellEnergy) > thr)
addNeighbour = true;
else
addNeighbour = false;
Expand All @@ -417,31 +426,31 @@ CaloTopoClusterFCCee::searchForNeighbours (
}
// If cell is hit.. but is assigned to another cluster
else if (itAllUsedCells != alreadyUsedCells.end() && itAllUsedCells->second != aClusterID && aAllowClusterMerge) {
uint32_t clusterIDToMerge = itAllUsedCells->second;
uint32_t clusterIDToMergeTo = itAllUsedCells->second;
if (msgLevel() <= MSG::VERBOSE){
verbose() << "This neighbour was found in cluster " << clusterIDToMerge
verbose() << "This neighbour was found in cluster " << clusterIDToMergeTo
<< ", cluster " << aClusterID
<< " will be merged!" << endmsg;
verbose() << "Assigning all cells ( "
<< protoClusters[aClusterID].size() << " ) to Cluster "
<< clusterIDToMerge << " with ( "
<< protoClusters[clusterIDToMerge].size()
<< clusterIDToMergeTo << " with ( "
<< protoClusters[clusterIDToMergeTo].size()
<< " ). " << endmsg;
}
// Fill all cells into cluster, and assigned cells to new cluster
alreadyUsedCells[neighbourID] = clusterIDToMerge;
alreadyUsedCells[neighbourID] = clusterIDToMergeTo;
for (const auto& cell : protoClusters[aClusterID]) {
alreadyUsedCells[cell.getCellID()] = clusterIDToMerge;
alreadyUsedCells[cell.getCellID()] = clusterIDToMergeTo;
// make sure that already assigned cells are not added
for (const auto& cellToMerge : protoClusters[clusterIDToMerge]) {
if (cellToMerge == cell) continue;
if (cellIdInColl(cell.getCellID(), protoClusters[clusterIDToMergeTo])) {
continue;
}
protoClusters[clusterIDToMerge].push_back(cell);
protoClusters[clusterIDToMergeTo].push_back(cell.clone());
}
protoClusters.erase(aClusterID);
// changed clusterId -> if more neighbours are found, correct assignment
verbose() << "Cluster Id changed to " << clusterIDToMerge << endmsg;
aClusterID = clusterIDToMerge;
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
Expand All @@ -463,7 +472,7 @@ void CaloTopoClusterFCCee::createCache(

// fill all noises and offsets values
for (const auto& cell : *aCells) {
int system = m_decoder->get(cell.getCellID(), m_index_system);
int system = m_decoder->get(cell.getCellID(), m_indexSystem);
if (system == 4) { //ECal barrel
int layer = m_decoder_ecal->get(cell.getCellID(), m_index_layer_ecal);
if (layers.find(layer) == layers.end()) {
Expand All @@ -488,3 +497,16 @@ void CaloTopoClusterFCCee::createCache(
m_min_noise[n.first] = *std::min_element(n.second.begin(), n.second.end());
}
}


inline bool CaloTopoClusterFCCee::cellIdInColl(
const uint64_t cellId,
const edm4hep::CalorimeterHitCollection& coll) {
for(const auto& cell: coll) {
if (cell.getCellID() == cellId) {
return true;
}
}

return false;
}
21 changes: 13 additions & 8 deletions RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,10 @@
#include "k4Interface/ICellPositionsTool.h"
#include "k4Interface/ITopoClusterInputTool.h"

#include <Gaudi/Property.h>
#include <cstdint>
#include <sys/types.h>
#include <vector>
#include <unordered_map>
#include <map>

class IGeoSvc;
Expand Down Expand Up @@ -115,7 +115,7 @@ class CaloTopoClusterFCCee : public GaudiAlgorithm {
/// Handle for neighbors tool
ToolHandle<ICaloReadNeighboursMap> m_neighboursTool{"TopoCaloNeighbours", this};

/// no segmentation used in HCal
/// No segmentation used in HCal
Gaudi::Property<bool> m_noSegmentationHCalUsed{this, "noSegmentationHCal", true, "HCal Barrel readout without DD4hep theta-module segmentation used."};
/// Seed threshold in sigma
Gaudi::Property<int> m_seedSigma{this, "seedSigma", 4, "number of sigma in noise threshold"};
Expand All @@ -126,22 +126,27 @@ class CaloTopoClusterFCCee : public GaudiAlgorithm {
/// Name of the electromagnetic calorimeter readout
Gaudi::Property<std::string> m_readoutName{this, "readoutName", "ECalBarrelModuleThetaMerged"};

/// 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<std::string> 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;

// Map of all calorimeter cells, key: cellID
std::map<uint64_t, const edm4hep::CalorimeterHit> m_allCellsMap;

// 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<double> m_min_offset;
std::vector<double> m_min_noise;

// Utility functions
void createCache(const edm4hep::CalorimeterHitCollection* aCells);
inline bool cellIdInColl(const uint64_t cellId,
const edm4hep::CalorimeterHitCollection& coll);

};
#endif /* RECFCCEECALORIMETER_CALOTOPOCLUSTERFCCEE_H */

0 comments on commit 71106d5

Please sign in to comment.