Skip to content

Commit

Permalink
Merge branch 'HEP-FCC:main' into gmarchio-main-20240724-positionedcells
Browse files Browse the repository at this point in the history
  • Loading branch information
giovannimarchiori authored Jul 25, 2024
2 parents 8f20000 + 05a64eb commit de7d18d
Show file tree
Hide file tree
Showing 19 changed files with 593 additions and 545 deletions.
76 changes: 44 additions & 32 deletions RecCalorimeter/src/components/ConstNoiseTool.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -17,29 +17,18 @@ ConstNoiseTool::ConstNoiseTool(const std::string& type, const std::string& name,
}

StatusCode ConstNoiseTool::initialize() {
// set for ECalBarrel system:5
// set for ECalEndcap system:6
// set for HCalEndcap system:7
// set for HCalBarrel system:8
// set for HCalExtBarrel system:9
// set for ECalFwd system:10
// set for HCalFwd system:11

m_systemNoiseConstMap.emplace(5, m_ECalBThreshold );
m_systemNoiseConstMap.emplace(6, m_ECalECThreshold );
m_systemNoiseConstMap.emplace(7, m_HCalECThreshold );
m_systemNoiseConstMap.emplace(8, m_HCalBThreshold );
m_systemNoiseConstMap.emplace(9, m_HCalBThreshold );
m_systemNoiseConstMap.emplace(10, m_ECalFwdThreshold );
m_systemNoiseConstMap.emplace(11, m_HCalFwdThreshold );

m_systemNoiseOffsetMap.emplace(5, 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. );

// check that sizes of m_detectors and m_detectorNoiseRMS and m_detectorNoiseOffset are the same
if (m_detectorsNoiseRMS.size() != m_detectors.size()) {
error() << "Sizes of the detectors vector and detectorsNoiseRMS vector do not match, exiting!" << endmsg;
return StatusCode::FAILURE;
}
if (m_setNoiseOffset) {
if (m_detectorsNoiseOffset.size() != m_detectors.size()) {
error() << "Sizes of the detectors vector and detectorsNoiseOffset vector do not match, exiting!" << endmsg;
return StatusCode::FAILURE;
}
}

// Get GeoSvc
m_geoSvc = service("GeoSvc");
Expand All @@ -49,6 +38,25 @@ StatusCode ConstNoiseTool::initialize() {
return StatusCode::FAILURE;
}

// loop over the detectors
for (size_t iDet=0; iDet<m_detectors.size(); iDet++) {
std::string detName = "DetID_" + m_detectors[iDet];
try {
int detID = m_geoSvc->getDetector()->constant<int>(detName);
m_systemNoiseRMSMap.emplace(detID, m_detectorsNoiseRMS[iDet]);
debug() << "Set noise RMS for detector " << detName << " (ID=" << detID << ") to: "
<< m_detectorsNoiseRMS[iDet] << endmsg;
if (m_setNoiseOffset) {
m_systemNoiseOffsetMap.emplace(detID, m_detectorsNoiseOffset[iDet]);
debug() << "Set noise offset for detector " << detName << " (ID=" << detID << ") to: "
<< m_detectorsNoiseOffset[iDet] << endmsg;
}
} catch (const std::exception& e) {
error() << "Error: detector with name " << detName << " not found, exception raised: " << e.what() << endmsg;
return StatusCode::FAILURE;
}
}

StatusCode sc = GaudiTool::initialize();
if (sc.isFailure()) return sc;

Expand All @@ -60,30 +68,34 @@ StatusCode ConstNoiseTool::finalize() {
return sc;
}

double ConstNoiseTool::getNoiseConstantPerCell(uint64_t aCellId) {
double ConstNoiseTool::getNoiseRMSPerCell(uint64_t aCellId) {

double Noise = 0.;
double noiseRMS = 0.;
// Get cells global coordinate "system"
dd4hep::DDSegmentation::CellID cID = aCellId;
unsigned cellSystem = m_decoder->get(cID, "system");
// cell noise in system
if (m_systemNoiseConstMap.count(cellSystem))
Noise = m_systemNoiseConstMap[cellSystem];
if (m_systemNoiseRMSMap.count(cellSystem))
noiseRMS = m_systemNoiseRMSMap[cellSystem];
else
warning() << "No noise constant set for subsystem " << cellSystem << "! Noise constant of cell set to 0. " << endmsg;
return Noise;
warning() << "No noise RMS set for subsystem " << cellSystem << "! Noise RMS of cell set to 0. " << endmsg;
return noiseRMS;
}

double ConstNoiseTool::getNoiseOffsetPerCell(uint64_t aCellId) {

double Noise = 0.;
if (!m_setNoiseOffset) {
return 0.;
}

double noiseOffset = 0.;
// Get cells global coordinate "system"
dd4hep::DDSegmentation::CellID cID = aCellId;
unsigned cellSystem = m_decoder->get(cID, "system");
// cell noise in system
if (m_systemNoiseOffsetMap.count(cellSystem))
Noise = m_systemNoiseOffsetMap[cellSystem];
noiseOffset = m_systemNoiseOffsetMap[cellSystem];
else
warning() << "No noise offset set for subsystem " << cellSystem << "! Noise offset of cell set to 0. " << endmsg;
return Noise;
return noiseOffset;
}
41 changes: 21 additions & 20 deletions RecCalorimeter/src/components/ConstNoiseTool.h
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@
#include "GaudiAlg/GaudiTool.h"
class IRndmGenSvc;

// k4geo
#include "detectorSegmentations/FCCSWGridPhiEta_k4geo.h"
// DD4HEP
#include "DDSegmentation/BitFieldCoder.h"

// k4FWCore
#include "k4Interface/INoiseConstTool.h"
Expand All @@ -16,12 +16,15 @@ class IGeoSvc;

/** @class ConstNoiseTool
*
* Tool to set noise constant for all cells in Calorimeters.
* Tool to set noise offset and RMS for all cells in Calorimeters.
* By default, set to rough estimates from ATLAS, can be changed in arguments for each Calo sub-system.
*
* @author Coralie Neubueser
* @date 2018-01
*
* @author Giovanni Marchiori
* @date 2024-07
*
*/

class ConstNoiseTool : public GaudiTool, virtual public INoiseConstTool {
Expand All @@ -34,34 +37,32 @@ class ConstNoiseTool : public GaudiTool, virtual public INoiseConstTool {
virtual StatusCode finalize() final;

/// Find the appropriate noise constant from the histogram
double getNoiseConstantPerCell(uint64_t aCellID);
double getNoiseRMSPerCell(uint64_t aCellID);
double getNoiseOffsetPerCell(uint64_t aCellID);

private:
std::map<uint,double> m_systemNoiseConstMap;
std::map<uint,double> m_systemNoiseRMSMap;
std::map<uint,double> m_systemNoiseOffsetMap;

/// Cell thresholds (1Sigma) in GeV:
/// List of subdetector names (they must match what is defined in DectDimension)
Gaudi::Property<std::vector<std::string>> m_detectors{this, "detectors", {"ECAL_Barrel", "ECAL_Endcap", "HCal_Endcap", "HCalBarrel", "HCalExtBarrel", "ECalFwd", "HCalFwd"}, "names of the calorimeters"};

/// Cell noise RMS (=1Sigma threshold) in GeV
/// default values estimates for the LAr and TileCal from ATLAS Barrel: ATL-LARG-PUB-2008_002
/// effective seed thresholds in topo-clustering of 7.5 and 11.5MeV in LAr and TileCal
/// EM Barrel Calorimeter
Gaudi::Property<double> m_ECalBThreshold{this, "ecalBarrelThreshold", 0.0075/4., "Cell threshold of EM Barrel Calorimeter in GeV"};
/// EM Endcap Calorimeter
Gaudi::Property<double> m_ECalECThreshold{this, "ecalEndcapThreshold", 0.0075/4., "Cell threshold of EM Endcap Calorimeter in GeV"};
/// EM Forward Calorimeter
Gaudi::Property<double> m_ECalFwdThreshold{this, "ecalFwdThreshold", 0.0075/4., "Cell threshold of EM Forward Calorimeter in GeV"};
/// Hadron Barrel Calorimeter
Gaudi::Property<double> m_HCalBThreshold{this, "hcalBarrelThreshold", 0.0115/4., "Cell threshold of Hadron Barrel Calorimeter in GeV"};
/// Hadron Endcap Calorimeter
Gaudi::Property<double> m_HCalECThreshold{this, "hcalEndcapThreshold", 0.0075/4., "Cell threshold of Hadron Endcap Calorimeter in GeV"};
/// Hadron Forward Calorimeter
Gaudi::Property<double> m_HCalFwdThreshold{this, "hcalFwdThreshold", 0.0075/4., "Cell threshold of Hadron Forward Calorimeter in GeV"};
Gaudi::Property<std::vector<double>> m_detectorsNoiseRMS{this, "detectorsNoiseRMS", {0.0075/4, 0.0075/4, 0.0115/4, 0.0115/4, 0.0115/4, 0.0075/4, 0.0075/4}, "Cell noise RMS in GeV"};

/// Set noise offset or not (false by default)
Gaudi::Property<bool> m_setNoiseOffset{this, "setNoiseOffset", false, "Set a noise offset per cell"};

/// Cell noise offset in GeV. 0 by default
Gaudi::Property<std::vector<double>> m_detectorsNoiseOffset{this, "detectorsNoiseOffset", {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, "Cell noise offset in GeV"};

/// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
// Decoder
dd4hep::DDSegmentation::BitFieldCoder* m_decoder = new dd4hep::DDSegmentation::BitFieldCoder("system:4");

/// Decoder
dd4hep::DDSegmentation::BitFieldCoder* m_decoder = new dd4hep::DDSegmentation::BitFieldCoder("system:4");
};

#endif /* RECCALORIMETER_CONSTNOISETOOL_H */
28 changes: 14 additions & 14 deletions RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -367,9 +367,9 @@ StatusCode CorrectECalBarrelSliWinCluster::execute() {
uint numCells = newCluster.hits_size();
double noise = 0;
if (m_constPileupNoise == 0) {
noise = getNoiseConstantPerCluster(newEta, numCells) * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
verbose() << " NUM CELLS = " << numCells << " cluster noise const = " << getNoiseConstantPerCluster(newEta, numCells)
<< " scaled to PU " << m_mu<< " = " << getNoiseConstantPerCluster(newEta, numCells)* std::sqrt(static_cast<int>(m_mu)) << endmsg;
noise = getNoiseRMSPerCluster(newEta, numCells) * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
verbose() << " NUM CELLS = " << numCells << " cluster noise const = " << getNoiseRMSPerCluster(newEta, numCells)
<< " scaled to PU " << m_mu<< " = " << getNoiseRMSPerCluster(newEta, numCells)* std::sqrt(static_cast<int>(m_mu)) << endmsg;
} else {
noise = m_constPileupNoise * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
}
Expand Down Expand Up @@ -456,36 +456,36 @@ StatusCode CorrectECalBarrelSliWinCluster::initNoiseFromFile() {
for (unsigned i = 0; i < 2; i++) {
pileupParamHistoName = m_pileupHistoName + std::to_string(i);
debug() << "Getting histogram with a name " << pileupParamHistoName << endmsg;
m_histoPileupConst.push_back(*dynamic_cast<TH1F*>(noiseFile->Get(pileupParamHistoName.c_str())));
if (m_histoPileupConst.at(i).GetNbinsX() < 1) {
m_histoPileupNoiseRMS.push_back(*dynamic_cast<TH1F*>(noiseFile->Get(pileupParamHistoName.c_str())));
if (m_histoPileupNoiseRMS.at(i).GetNbinsX() < 1) {
error() << "Histogram " << pileupParamHistoName
<< " has 0 bins! check the file with noise and the name of the histogram!" << endmsg;
return StatusCode::FAILURE;
}
}

// Check if we have same number of histograms (all layers) for pileup and electronics noise
if (m_histoPileupConst.size() == 0) {
if (m_histoPileupNoiseRMS.size() == 0) {
error() << "No histograms with noise found!!!!" << endmsg;
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
}

double CorrectECalBarrelSliWinCluster::getNoiseConstantPerCluster(double aEta, uint aNumCells) {
double CorrectECalBarrelSliWinCluster::getNoiseRMSPerCluster(double aEta, uint aNumCells) {
double param0 = 0.;
double param1 = 0.;

// All histograms have same binning, all bins with same size
// Using the histogram of the first parameter to get the bin size
unsigned index = 0;
if (m_histoPileupConst.size() != 0) {
int Nbins = m_histoPileupConst.at(index).GetNbinsX();
if (m_histoPileupNoiseRMS.size() != 0) {
int Nbins = m_histoPileupNoiseRMS.at(index).GetNbinsX();
double deltaEtaBin =
(m_histoPileupConst.at(index).GetBinLowEdge(Nbins) + m_histoPileupConst.at(index).GetBinWidth(Nbins) -
m_histoPileupConst.at(index).GetBinLowEdge(1)) /
(m_histoPileupNoiseRMS.at(index).GetBinLowEdge(Nbins) + m_histoPileupNoiseRMS.at(index).GetBinWidth(Nbins) -
m_histoPileupNoiseRMS.at(index).GetBinLowEdge(1)) /
Nbins;
double etaFirtsBin = m_histoPileupConst.at(index).GetBinLowEdge(1);
double etaFirtsBin = m_histoPileupNoiseRMS.at(index).GetBinLowEdge(1);
// find the eta bin for the cell
int ibin = floor((fabs(aEta) - etaFirtsBin) / deltaEtaBin) + 1;
verbose() << "Current eta = " << aEta << " bin = " << ibin << endmsg;
Expand All @@ -494,8 +494,8 @@ double CorrectECalBarrelSliWinCluster::getNoiseConstantPerCluster(double aEta, u
<< endmsg;
ibin = Nbins;
}
param0 = m_histoPileupConst.at(0).GetBinContent(ibin);
param1 = m_histoPileupConst.at(1).GetBinContent(ibin);
param0 = m_histoPileupNoiseRMS.at(0).GetBinContent(ibin);
param1 = m_histoPileupNoiseRMS.at(1).GetBinContent(ibin);
verbose() << "p0 = " << param0 << " param1 = " << param1 << endmsg;
} else {
debug() << "No histograms with noise constants!!!!! " << endmsg;
Expand Down
10 changes: 5 additions & 5 deletions RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.h
Original file line number Diff line number Diff line change
Expand Up @@ -81,12 +81,12 @@ class CorrectECalBarrelSliWinCluster : public GaudiAlgorithm {
* @return Status code if retriving histograms was successful
*/
StatusCode initNoiseFromFile();
/** Find the appropriate noise constant from the histogram
/** Find the appropriate noise RMS from the histogram
* @param[in] aEta Pseudorapidity value of the centre of the cluster
* @param[in] aNumCells Number of cells in a cluster
* @return Width of the Gaussian distribution of noise per cluster
* @return Width of the Gaussian distribution of noise per cluster
*/
double getNoiseConstantPerCluster(double aEta, uint numCells);
double getNoiseRMSPerCluster(double aEta, uint numCells);
/// Handle for clusters (input collection)
DataHandle<edm4hep::ClusterCollection> m_inClusters{"clusters", Gaudi::DataHandle::Reader, this};
/// Handle for corrected clusters (output collection)
Expand Down Expand Up @@ -171,8 +171,8 @@ class CorrectECalBarrelSliWinCluster : public GaudiAlgorithm {
"Name of pileup histogram"};
/// Pile-up scenario for pileup noise scaling (with sqrt(mu))
Gaudi::Property<int> m_mu{this, "mu", 1000, "Pileup scenario"};
/// Histograms with pileup constants (index in array - parameter number: 0 or 1)
std::vector<TH1F> m_histoPileupConst;
/// Histograms with pileup RMS (index in array - parameter number: 0 or 1)
std::vector<TH1F> m_histoPileupNoiseRMS;
/// Values of eta corresponding to the upstream correction parameters
Gaudi::Property<std::vector<double>> m_etaValues{this, "etaValues", {0}, "Eta values"};
/// Borders of the eta bins for the upstream correction (middle between eta values)
Expand Down
8 changes: 4 additions & 4 deletions RecCalorimeter/src/components/MassInv.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -545,9 +545,9 @@ StatusCode MassInv::execute() {
uint numCells = newCluster.hits_size();
double noise = 0;
if (m_constPileupNoise == 0) {
noise = getNoiseConstantPerCluster(newEta, numCells) * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
verbose() << " NUM CELLS = " << numCells << " cluster noise const = " << getNoiseConstantPerCluster(newEta, numCells)
<< " scaled to PU " << m_mu<< " = " << getNoiseConstantPerCluster(newEta, numCells)* std::sqrt(static_cast<int>(m_mu)) << endmsg;
noise = getNoiseRMSPerCluster(newEta, numCells) * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
verbose() << " NUM CELLS = " << numCells << " cluster noise RMS = " << getNoiseRMSPerCluster(newEta, numCells)
<< " scaled to PU " << m_mu<< " = " << noise << endmsg;
} else {
noise = m_constPileupNoise * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
}
Expand Down Expand Up @@ -849,7 +849,7 @@ StatusCode MassInv::initNoiseFromFile() {
return StatusCode::SUCCESS;
}

double MassInv::getNoiseConstantPerCluster(double aEta, uint aNumCells) {
double MassInv::getNoiseRMSPerCluster(double aEta, uint aNumCells) {
double param0 = 0.;
double param1 = 0.;

Expand Down
2 changes: 1 addition & 1 deletion RecCalorimeter/src/components/MassInv.h
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ class MassInv : public GaudiAlgorithm {
* @param[in] aNumCells Number of cells in a cluster
* @return Width of the Gaussian distribution of noise per cluster
*/
double getNoiseConstantPerCluster(double aEta, uint numCells);
double getNoiseRMSPerCluster(double aEta, uint numCells);
/// Handle for clusters (input collection)
DataHandle<edm4hep::ClusterCollection> m_inClusters{"clusters", Gaudi::DataHandle::Reader, this};
/// Handle for corrected clusters (output collection)
Expand Down
7 changes: 4 additions & 3 deletions RecCalorimeter/src/components/NoiseCaloCellsFlatTool.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -27,19 +27,20 @@ StatusCode NoiseCaloCellsFlatTool::initialize() {
}
}

info() << "Sigma of the cell noise: " << m_cellNoise * 1.e3 << " MeV" << endmsg;
info() << "RMS of the cell noise: " << m_cellNoiseRMS * 1.e3 << " MeV" << endmsg;
info() << "Offset of the cell noise: " << m_cellNoiseOffset * 1.e3 << " MeV" << endmsg;
info() << "Filter noise threshold: " << m_filterThreshold << "*sigma" << endmsg;
return StatusCode::SUCCESS;
}

void NoiseCaloCellsFlatTool::addRandomCellNoise(std::unordered_map<uint64_t, double>& aCells) {
std::for_each(aCells.begin(), aCells.end(),
[this](std::pair<const uint64_t, double>& p) { p.second += (m_gauss.shoot() * m_cellNoise); });
[this](std::pair<const uint64_t, double>& p) { p.second += (m_cellNoiseOffset + (m_gauss.shoot() * m_cellNoiseRMS)); });
}

void NoiseCaloCellsFlatTool::filterCellNoise(std::unordered_map<uint64_t, double>& aCells) {
// Erase a cell if it has energy below a threshold
double threshold = m_filterThreshold * m_cellNoise;
double threshold = m_cellNoiseOffset + m_filterThreshold * m_cellNoiseRMS;
auto it = aCells.begin();
while ((it = std::find_if(it, aCells.end(), [&threshold](std::pair<const uint64_t, double>& p) {
return bool(p.second < threshold);
Expand Down
16 changes: 10 additions & 6 deletions RecCalorimeter/src/components/NoiseCaloCellsFlatTool.h
Original file line number Diff line number Diff line change
Expand Up @@ -13,11 +13,13 @@
*
* Very simple tool for calorimeter noise using a single noise value for all cells
* createRandomCellNoise: Create random CaloHits (gaussian distribution) for the vector of cells
* filterCellNoise: remove cells with energy bellow threshold*sigma from the vector of cells
* filterCellNoise: remove cells with energy below threshold*sigma from the vector of cells
*
* @author Jana Faltova
* @date 2016-09
*
* @author Giovanni Marchiori
* @date 2024-07
*/

class NoiseCaloCellsFlatTool : public GaudiTool, virtual public INoiseCaloCellsTool {
Expand All @@ -31,17 +33,19 @@ class NoiseCaloCellsFlatTool : public GaudiTool, virtual public INoiseCaloCellsT
* Vector of cells must contain all cells in the calorimeter with their cellIDs.
*/
virtual void addRandomCellNoise(std::unordered_map<uint64_t, double>& aCells) final;
/** @brief Remove cells with energy bellow threshold*sigma from the vector of cells
/** @brief Remove cells with energy below threshold*sigma from the vector of cells
*/
virtual void filterCellNoise(std::unordered_map<uint64_t, double>& aCells) final;

private:
/// Sigma of noise -- uniform noise per cell in GeV
Gaudi::Property<double> m_cellNoise{this, "cellNoise", 0.003, "uniform noise per cell in GeV"};
/// Energy threshold (Ecell < filterThreshold*m_cellNoise removed)
/// RMS of noise -- uniform RMS per cell in GeV
Gaudi::Property<double> m_cellNoiseRMS{this, "cellNoiseRMS", 0.003, "uniform noise RMS per cell in GeV"};
/// Offset of noise -- uniform offset per cell in GeV
Gaudi::Property<double> m_cellNoiseOffset{this, "cellNoiseOffset", 0.0, "uniform noise offset per cell in GeV"};
/// Energy threshold (Ecell < m_cellNoiseOffset + filterThreshold*m_cellNoiseRMS removed)
Gaudi::Property<double> m_filterThreshold{
this, "filterNoiseThreshold", 3,
"remove cells with energy bellow filterThreshold (threshold is multiplied by a cell noise sigma)"};
"remove cells with energy below offset + threshold * noise RMS"};
/// Random Number Service
IRndmGenSvc* m_randSvc;
/// Gaussian random number generator used for smearing with a constant resolution (m_sigma)
Expand Down
Loading

0 comments on commit de7d18d

Please sign in to comment.