Skip to content

Commit

Permalink
make topoclustering work with new module-theta merged readout
Browse files Browse the repository at this point in the history
  • Loading branch information
@giovannimarchiori
giovannimarchiori committed Sep 27, 2023
1 parent e1f62eb commit e633941
Show file tree
Hide file tree
Showing 9 changed files with 977 additions and 35 deletions.
54 changes: 44 additions & 10 deletions RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ DECLARE_COMPONENT(ReadNoiseFromFileTool)
ReadNoiseFromFileTool::ReadNoiseFromFileTool(const std::string& type, const std::string& name, const IInterface* parent)
: GaudiTool(type, name, parent) {
declareInterface<INoiseConstTool>(this);
declareProperty("cellPositionsTool", m_cellPositionsTool, "Handle for tool to retrieve cell positions");
}

StatusCode ReadNoiseFromFileTool::initialize() {
Expand All @@ -29,10 +30,24 @@ StatusCode ReadNoiseFromFileTool::initialize() {
<< "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
return StatusCode::FAILURE;
}
m_segmentation = dynamic_cast<dd4hep::DDSegmentation::FCCSWGridPhiEta*>(m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation());
if (m_segmentation == nullptr) {
error() << "There is no phi-eta segmentation!!!!" << endmsg;
return StatusCode::FAILURE;

// Check if cell position tool available if m_useSeg==false; if tool not
// available, try using segmentation instead
if (!m_cellPositionsTool.retrieve() and !m_useSeg) {
info() << "Unable to retrieve cell positions tool, try eta-phi segmentation." << endmsg;
m_useSeg = true;
}

// Get PhiEta segmentation
if (m_useSeg) {
m_segmentation = dynamic_cast<dd4hep::DDSegmentation::FCCSWGridPhiEta*>(
m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation());
if (m_segmentation == nullptr) {
error() << "There is no phi-eta segmentation." << endmsg;
return StatusCode::FAILURE;
}
else
info() << "Found phi-eta segmentation." << endmsg;
}

// open and check file, read the histograms with noise constants
Expand All @@ -43,6 +58,10 @@ StatusCode ReadNoiseFromFileTool::initialize() {

// Take readout bitfield decoder from GeoSvc
m_decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder();
if (m_decoder == nullptr) {
error() << "Cannot create decore for readout " << m_readoutName << endmsg;
return StatusCode::FAILURE;
}

StatusCode sc = GaudiTool::initialize();
if (sc.isFailure()) return sc;
Expand Down Expand Up @@ -134,21 +153,28 @@ double ReadNoiseFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) {

// Get cell coordinates: eta and radial layer
dd4hep::DDSegmentation::CellID cID = aCellId;
double cellEta = m_segmentation->eta(cID);

double cellEta;
if (m_useSeg)
cellEta = m_segmentation->eta(aCellId);
else
cellEta = m_cellPositionsTool->xyzPosition(aCellId).Eta();
unsigned cellLayer = m_decoder->get(cID, m_activeFieldName);

// All histograms have same binning, all bins with same size
// Using the histogram in the first layer to get the bin size
unsigned index = 0;
if (m_histoElecNoiseConst.size() != 0) {
int Nbins = m_histoElecNoiseConst.at(index).GetNbinsX();
// GM: basically the same as FindBin ...
/*
double deltaEtaBin =
(m_histoElecNoiseConst.at(index).GetBinLowEdge(Nbins) + m_histoElecNoiseConst.at(index).GetBinWidth(Nbins) -
m_histoElecNoiseConst.at(index).GetBinLowEdge(1)) /
Nbins;
// find the eta bin for the cell
int ibin = floor(fabs(cellEta) / deltaEtaBin) + 1;
*/
int ibin = m_histoElecNoiseConst.at(index).FindBin(fabs(cellEta));
if (ibin > Nbins) {
error() << "eta outside range of the histograms! Cell eta: " << cellEta << " Nbins in histogram: " << Nbins
<< endmsg;
Expand Down Expand Up @@ -189,20 +215,28 @@ double ReadNoiseFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) {

// Get cell coordinates: eta and radial layer
dd4hep::DDSegmentation::CellID cID = aCellId;
double cellEta = m_segmentation->eta(cID);
double cellEta;
if (m_useSeg)
cellEta = m_segmentation->eta(aCellId);
else
cellEta = m_cellPositionsTool->xyzPosition(aCellId).Eta();
unsigned cellLayer = m_decoder->get(cID, m_activeFieldName);

// All histograms have same binning, all bins with same size
// Using the histogram in the first layer to get the bin size
unsigned index = 0;
if (m_histoElecNoiseOffset.size() != 0) {
int Nbins = m_histoElecNoiseOffset.at(index).GetNbinsX();
double deltaEtaBin =
// find the eta bin for the cell
// GM: basically the same as FindBin ...
/*
double deltaEtaBin =
(m_histoElecNoiseOffset.at(index).GetBinLowEdge(Nbins) + m_histoElecNoiseOffset.at(index).GetBinWidth(Nbins) -
m_histoElecNoiseOffset.at(index).GetBinLowEdge(1)) /
Nbins;
// find the eta bin for the cell
int ibin = floor(fabs(cellEta) / deltaEtaBin) + 1;
int ibin = floor(fabs(cellEta) / deltaEtaBin) + 1;
*/
int ibin = m_histoElecNoiseOffset.at(index).FindBin(fabs(cellEta));
if (ibin > Nbins) {
error() << "eta outside range of the histograms! Cell eta: " << cellEta << " Nbins in histogram: " << Nbins
<< endmsg;
Expand Down
17 changes: 8 additions & 9 deletions RecCalorimeter/src/components/ReadNoiseFromFileTool.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,12 @@

// FCCSW
#include "k4FWCore/DataHandle.h"
#include "DetSegmentation/FCCSWGridPhiEta.h"
#include "k4Interface/ICalorimeterTool.h"
#include "k4Interface/INoiseConstTool.h"
#include "k4Interface/ICellPositionsTool.h"

#include "DetSegmentation/FCCSWGridPhiEta.h"

class IGeoSvc;

// Root
Expand Down Expand Up @@ -43,15 +44,18 @@ class ReadNoiseFromFileTool : public GaudiTool, virtual public INoiseConstTool {
double getNoiseOffsetPerCell(uint64_t aCellID);

private:
/// Handle for tool to get cell positions
ToolHandle<ICellPositionsTool> m_cellPositionsTool;
/// use segmentation (eta-phi only so far!) in case no cell position tool is defined
Gaudi::Property<bool> m_useSeg{this, "useSegmentation", true, "Specify if segmentation is used to determine cell position."};
/// Add pileup contribution to the electronics noise? (only if read from file)
Gaudi::Property<bool> m_addPileup{this, "addPileup", true,
"Add pileup contribution to the electronics noise? (only if read from file)"};
/// Noise offset, if false, mean is set to 0
Gaudi::Property<bool> m_setNoiseOffset{this, "setNoiseOffset", true, "Set a noise offset per cell"};

/// Name of the file with noise constants
Gaudi::Property<std::string> m_noiseFileName{this, "noiseFileName", "", "Name of the file with noise constants"};
/// Name of the detector readout
/// Name of the detector readout (needed to get the decoder to extract e.g. layer information)
Gaudi::Property<std::string> m_readoutName{this, "readoutName", "ECalHitsPhiEta", "Name of the detector readout"};
/// Name of active layers for sampling calorimeter
Gaudi::Property<std::string> m_activeFieldName{this, "activeFieldName", "active_layer",
Expand All @@ -66,28 +70,23 @@ class ReadNoiseFromFileTool : public GaudiTool, virtual public INoiseConstTool {
"Name of electronics noise offset histogram"};
/// Name of pileup offset histogram
Gaudi::Property<std::string> m_pileupOffsetHistoName{this, "pileupOffsetHistoName", "h_pileup_layer", "Name of pileup offset histogram"};

/// Number of radial layers
Gaudi::Property<uint> m_numRadialLayers{this, "numRadialLayers", 3, "Number of radial layers"};

/// Factor to apply to the noise values to get them in GeV if e.g. they were produced in MeV
Gaudi::Property<float> m_scaleFactor{this, "scaleFactor", 1, "Factor to apply to the noise values"};

/// Histograms with pileup constants (index in array - radial layer)
std::vector<TH1F> m_histoPileupConst;
/// Histograms with electronics noise constants (index in array - radial layer)
std::vector<TH1F> m_histoElecNoiseConst;

/// Histograms with pileup offset (index in array - radial layer)
std::vector<TH1F> m_histoPileupOffset;
/// Histograms with electronics noise offset (index in array - radial layer)
std::vector<TH1F> m_histoElecNoiseOffset;

/// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
/// PhiEta segmentation
dd4hep::DDSegmentation::FCCSWGridPhiEta* m_segmentation;
// Decoder
// Decoder for ECal layers
dd4hep::DDSegmentation::BitFieldCoder* m_decoder;

};
Expand Down
1 change: 1 addition & 0 deletions RecFCCeeCalorimeter/CMakeLists.txt
View file
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
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
4 changes: 2 additions & 2 deletions RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -126,15 +126,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

0 comments on commit e633941

Please sign in to comment.