Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

make topoclustering work with new module-theta merged readout #51

Merged
merged 6 commits into from
Oct 24, 2023
Merged

make topoclustering work with new module-theta merged readout #51

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
e633941
make topoclustering work with new module-theta merged readout
giovannimarchiori Sep 27, 2023
b491630
fix for theta readout
giovannimarchiori Sep 27, 2023
68792c4
move ReadNoiseFromFileTool for theta-based readout to new tool in FCC…
giovannimarchiori Sep 29, 2023
72b9814
make inclusion of diagonal neighbours configurable
giovannimarchiori Oct 9, 2023
a848093
rename value of a property
giovannimarchiori Oct 9, 2023
66d6e75
check that tools are retrieved successfully
giovannimarchiori Oct 24, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions RecFCCeeCalorimeter/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ gaudi_add_module(k4RecFCCeeCalorimeterPlugins
DD4hep::DDCore
EDM4HEP::edm4hep
FCCDetectors::DetSegmentation
FCCDetectors::DetCommon
DD4hep::DDG4
ROOT::Core
ROOT::Hist
Expand Down
43 changes: 29 additions & 14 deletions RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -117,8 +117,15 @@ StatusCode CaloTopoClusterFCCee::execute() {
});

std::map<uint, std::vector<std::pair<uint64_t, int>>> preClusterCollection;
CaloTopoClusterFCCee::buildingProtoCluster(m_neighbourSigma, m_lastNeighbourSigma, firstSeeds, m_allCells,
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;
}

// Build Clusters in edm
debug() << "Building " << preClusterCollection.size() << " cluster." << endmsg;
double checkTotEnergy = 0.;
Expand All @@ -132,10 +139,10 @@ StatusCode CaloTopoClusterFCCee::execute() {
double energy = 0.;
double deltaR = 0.;
std::vector<double> posPhi (i.second.size());
std::vector<double> posEta (i.second.size());
std::vector<double> posTheta (i.second.size());
std::vector<double> vecEnergy (i.second.size());
double sumPhi = 0.;
double sumEta = 0.;
double sumTheta = 0.;
std::map<int,int> system;

for (auto pair : i.second) {
Expand Down Expand Up @@ -178,10 +185,10 @@ StatusCode CaloTopoClusterFCCee::execute() {
posY += posCell.Y() * newCell.getEnergy();
posZ += posCell.Z() * newCell.getEnergy();
posPhi.push_back(posCell.Phi());
posEta.push_back(posCell.Eta());
posTheta.push_back(posCell.Theta());
vecEnergy.push_back(newCell.getEnergy());
sumPhi += posCell.Phi() * newCell.getEnergy();
sumEta += posCell.Eta() * newCell.getEnergy();
sumTheta += posCell.Theta() * newCell.getEnergy();

cluster.addToHits(newCell);
auto er = m_allCells.erase(cID);
Expand All @@ -193,10 +200,10 @@ StatusCode CaloTopoClusterFCCee::execute() {
cluster.setPosition(edm4hep::Vector3f(posX / energy, posY / energy, posZ / energy));
// store deltaR of cluster in time for the moment..
sumPhi = sumPhi / energy;
sumEta = sumEta / energy;
sumTheta = sumTheta / energy;
int counter = 0;
for (auto entryEta : posEta){
deltaR += sqrt(pow(entryEta-sumEta,2) + pow(posEta[counter]-sumPhi,2)) * vecEnergy[counter];
for (auto entryTheta : posTheta){
deltaR += sqrt(pow(entryTheta-sumTheta,2) + pow(posPhi[counter]-sumPhi,2)) * vecEnergy[counter];
counter++;
}
cluster.addToShapeParameters(deltaR / energy);
Expand All @@ -208,7 +215,7 @@ StatusCode CaloTopoClusterFCCee::execute() {
clusterWithMixedCells++;

posPhi.clear();
posEta.clear();
posTheta.clear();
vecEnergy.clear();

}
Expand All @@ -225,12 +232,19 @@ void CaloTopoClusterFCCee::findingSeeds(const std::unordered_map<uint64_t, doubl
std::vector<std::pair<uint64_t, double>>& aSeeds) {
for (const auto& cell : aCells) {
// 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);

double min_threshold = m_min_offset[layer] + m_min_noise[layer] * aNumSigma;

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() << "min_threshold = " << min_threshold << endmsg;
debug() << "abs(cell.second) = " << abs(cell.second) << 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
continue;
Expand All @@ -240,9 +254,10 @@ void CaloTopoClusterFCCee::findingSeeds(const std::unordered_map<uint64_t, doubl
// 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){
verbose() << "noise offset = " << m_noiseTool->noiseOffset(cell.first) << "GeV " << endmsg;
verbose() << "noise rms = " << m_noiseTool->noiseRMS(cell.first) << "GeV " << endmsg;
verbose() << "seed threshold = " << threshold << "GeV " << endmsg;
debug() << "noise offset = " << m_noiseTool->noiseOffset(cell.first) << "GeV " << endmsg;
debug() << "noise rms = " << m_noiseTool->noiseRMS(cell.first) << "GeV " << endmsg;
debug() << "seed threshold = " << threshold << "GeV " << endmsg;
debug() << "======================================" << endmsg;
}
if (abs(cell.second) > threshold) {
aSeeds.emplace_back(cell.first, cell.second);
Expand Down
7 changes: 3 additions & 4 deletions RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,7 @@ class Segmentation;
}
}

/** @class CaloTopoClusterFCCeeAlgorithm Reconstruction/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h
* CombinedCaloTopoCluster.h
/** @class CaloTopoClusterFCCee k4RecCalorimeter/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h
*
* Algorithm building the topological clusters for the energy reconstruction, following ATLAS note
* ATL-LARG-PUB-2008-002.
Expand Down Expand Up @@ -126,15 +125,15 @@ class CaloTopoClusterFCCee : public GaudiAlgorithm {
//ToolHandle<ICellPositionsTool> m_cellPositionsHFwdTool{"CellPositionsCaloDiscsTool", this};

/// no segmentation used in HCal
Gaudi::Property<bool> m_noSegmentationHCalUsed{this, "noSegmentationHCal", true, "HCal Barrel readout without DD4hep eta-phi segmentation used."};
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"};
/// Neighbour threshold in sigma
Gaudi::Property<int> m_neighbourSigma{this, "neighbourSigma", 2, "number of sigma in noise threshold"};
/// Last neighbour threshold in sigma
Gaudi::Property<int> m_lastNeighbourSigma{this, "lastNeighbourSigma", 0, "number of sigma in noise threshold"};
/// Name of the electromagnetic calorimeter readout
Gaudi::Property<std::string> m_readoutName{this, "readoutName", "ECalBarrelPhiEta"};
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");
Expand Down
Loading
Loading