diff --git a/RecCalorimeter/src/components/ConstNoiseTool.cpp b/RecCalorimeter/src/components/ConstNoiseTool.cpp index dc007dde..f446d4cb 100644 --- a/RecCalorimeter/src/components/ConstNoiseTool.cpp +++ b/RecCalorimeter/src/components/ConstNoiseTool.cpp @@ -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"); @@ -49,6 +38,25 @@ StatusCode ConstNoiseTool::initialize() { return StatusCode::FAILURE; } + // loop over the detectors + for (size_t iDet=0; iDetgetDetector()->constant(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; @@ -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; } diff --git a/RecCalorimeter/src/components/ConstNoiseTool.h b/RecCalorimeter/src/components/ConstNoiseTool.h index bc62195f..b4f2c798 100644 --- a/RecCalorimeter/src/components/ConstNoiseTool.h +++ b/RecCalorimeter/src/components/ConstNoiseTool.h @@ -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" @@ -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 { @@ -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 m_systemNoiseConstMap; + std::map m_systemNoiseRMSMap; std::map m_systemNoiseOffsetMap; - /// Cell thresholds (1Sigma) in GeV: + /// List of subdetector names (they must match what is defined in DectDimension) + Gaudi::Property> 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 m_ECalBThreshold{this, "ecalBarrelThreshold", 0.0075/4., "Cell threshold of EM Barrel Calorimeter in GeV"}; - /// EM Endcap Calorimeter - Gaudi::Property m_ECalECThreshold{this, "ecalEndcapThreshold", 0.0075/4., "Cell threshold of EM Endcap Calorimeter in GeV"}; - /// EM Forward Calorimeter - Gaudi::Property m_ECalFwdThreshold{this, "ecalFwdThreshold", 0.0075/4., "Cell threshold of EM Forward Calorimeter in GeV"}; - /// Hadron Barrel Calorimeter - Gaudi::Property m_HCalBThreshold{this, "hcalBarrelThreshold", 0.0115/4., "Cell threshold of Hadron Barrel Calorimeter in GeV"}; - /// Hadron Endcap Calorimeter - Gaudi::Property m_HCalECThreshold{this, "hcalEndcapThreshold", 0.0075/4., "Cell threshold of Hadron Endcap Calorimeter in GeV"}; - /// Hadron Forward Calorimeter - Gaudi::Property m_HCalFwdThreshold{this, "hcalFwdThreshold", 0.0075/4., "Cell threshold of Hadron Forward Calorimeter in GeV"}; + Gaudi::Property> 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 m_setNoiseOffset{this, "setNoiseOffset", false, "Set a noise offset per cell"}; + + /// Cell noise offset in GeV. 0 by default + Gaudi::Property> 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 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 */ diff --git a/RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.cpp b/RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.cpp index ee0f19e3..502a58ab 100644 --- a/RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.cpp +++ b/RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.cpp @@ -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(m_mu)); - verbose() << " NUM CELLS = " << numCells << " cluster noise const = " << getNoiseConstantPerCluster(newEta, numCells) - << " scaled to PU " << m_mu<< " = " << getNoiseConstantPerCluster(newEta, numCells)* std::sqrt(static_cast(m_mu)) << endmsg; + noise = getNoiseRMSPerCluster(newEta, numCells) * m_gauss.shoot() * std::sqrt(static_cast(m_mu)); + verbose() << " NUM CELLS = " << numCells << " cluster noise const = " << getNoiseRMSPerCluster(newEta, numCells) + << " scaled to PU " << m_mu<< " = " << getNoiseRMSPerCluster(newEta, numCells)* std::sqrt(static_cast(m_mu)) << endmsg; } else { noise = m_constPileupNoise * m_gauss.shoot() * std::sqrt(static_cast(m_mu)); } @@ -456,8 +456,8 @@ 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(noiseFile->Get(pileupParamHistoName.c_str()))); - if (m_histoPileupConst.at(i).GetNbinsX() < 1) { + m_histoPileupNoiseRMS.push_back(*dynamic_cast(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; @@ -465,27 +465,27 @@ StatusCode CorrectECalBarrelSliWinCluster::initNoiseFromFile() { } // 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; @@ -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; diff --git a/RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.h b/RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.h index f26401fc..e0174de4 100644 --- a/RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.h +++ b/RecCalorimeter/src/components/CorrectECalBarrelSliWinCluster.h @@ -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 m_inClusters{"clusters", Gaudi::DataHandle::Reader, this}; /// Handle for corrected clusters (output collection) @@ -171,8 +171,8 @@ class CorrectECalBarrelSliWinCluster : public GaudiAlgorithm { "Name of pileup histogram"}; /// Pile-up scenario for pileup noise scaling (with sqrt(mu)) Gaudi::Property m_mu{this, "mu", 1000, "Pileup scenario"}; - /// Histograms with pileup constants (index in array - parameter number: 0 or 1) - std::vector m_histoPileupConst; + /// Histograms with pileup RMS (index in array - parameter number: 0 or 1) + std::vector m_histoPileupNoiseRMS; /// Values of eta corresponding to the upstream correction parameters Gaudi::Property> m_etaValues{this, "etaValues", {0}, "Eta values"}; /// Borders of the eta bins for the upstream correction (middle between eta values) diff --git a/RecCalorimeter/src/components/MassInv.cpp b/RecCalorimeter/src/components/MassInv.cpp index d38c14c2..707684df 100644 --- a/RecCalorimeter/src/components/MassInv.cpp +++ b/RecCalorimeter/src/components/MassInv.cpp @@ -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(m_mu)); - verbose() << " NUM CELLS = " << numCells << " cluster noise const = " << getNoiseConstantPerCluster(newEta, numCells) - << " scaled to PU " << m_mu<< " = " << getNoiseConstantPerCluster(newEta, numCells)* std::sqrt(static_cast(m_mu)) << endmsg; + noise = getNoiseRMSPerCluster(newEta, numCells) * m_gauss.shoot() * std::sqrt(static_cast(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(m_mu)); } @@ -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.; diff --git a/RecCalorimeter/src/components/MassInv.h b/RecCalorimeter/src/components/MassInv.h index 833b6763..28ae4df9 100644 --- a/RecCalorimeter/src/components/MassInv.h +++ b/RecCalorimeter/src/components/MassInv.h @@ -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 m_inClusters{"clusters", Gaudi::DataHandle::Reader, this}; /// Handle for corrected clusters (output collection) diff --git a/RecCalorimeter/src/components/NoiseCaloCellsFlatTool.cpp b/RecCalorimeter/src/components/NoiseCaloCellsFlatTool.cpp index ea36f23a..0e779adc 100644 --- a/RecCalorimeter/src/components/NoiseCaloCellsFlatTool.cpp +++ b/RecCalorimeter/src/components/NoiseCaloCellsFlatTool.cpp @@ -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& aCells) { std::for_each(aCells.begin(), aCells.end(), - [this](std::pair& p) { p.second += (m_gauss.shoot() * m_cellNoise); }); + [this](std::pair& p) { p.second += (m_cellNoiseOffset + (m_gauss.shoot() * m_cellNoiseRMS)); }); } void NoiseCaloCellsFlatTool::filterCellNoise(std::unordered_map& 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& p) { return bool(p.second < threshold); diff --git a/RecCalorimeter/src/components/NoiseCaloCellsFlatTool.h b/RecCalorimeter/src/components/NoiseCaloCellsFlatTool.h index aa5ee8f3..2e170efb 100644 --- a/RecCalorimeter/src/components/NoiseCaloCellsFlatTool.h +++ b/RecCalorimeter/src/components/NoiseCaloCellsFlatTool.h @@ -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 { @@ -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& 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& aCells) final; private: - /// Sigma of noise -- uniform noise per cell in GeV - Gaudi::Property 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 m_cellNoiseRMS{this, "cellNoiseRMS", 0.003, "uniform noise RMS per cell in GeV"}; + /// Offset of noise -- uniform offset per cell in GeV + Gaudi::Property m_cellNoiseOffset{this, "cellNoiseOffset", 0.0, "uniform noise offset per cell in GeV"}; + /// Energy threshold (Ecell < m_cellNoiseOffset + filterThreshold*m_cellNoiseRMS removed) Gaudi::Property 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) diff --git a/RecCalorimeter/src/components/NoiseCaloCellsFromFileTool.cpp b/RecCalorimeter/src/components/NoiseCaloCellsFromFileTool.cpp index cf23a7fa..424a30bb 100644 --- a/RecCalorimeter/src/components/NoiseCaloCellsFromFileTool.cpp +++ b/RecCalorimeter/src/components/NoiseCaloCellsFromFileTool.cpp @@ -97,7 +97,7 @@ StatusCode NoiseCaloCellsFromFileTool::initialize() { void NoiseCaloCellsFromFileTool::addRandomCellNoise(std::unordered_map& aCells) { std::for_each(aCells.begin(), aCells.end(), [this](std::pair& p) { - p.second += (getNoiseConstantPerCell(p.first) * m_gauss.shoot()); + p.second += (getNoiseRMSPerCell(p.first) * m_gauss.shoot()); }); } @@ -105,7 +105,7 @@ void NoiseCaloCellsFromFileTool::filterCellNoise(std::unordered_map& p) { - return m_useAbsInFilter ? bool(std::abs(p.second) < m_filterThreshold * getNoiseConstantPerCell(p.first)) : bool(p.second < m_filterThreshold * getNoiseConstantPerCell(p.first)); + return m_useAbsInFilter ? bool(std::abs(p.second) < m_filterThreshold * getNoiseRMSPerCell(p.first)) : bool(p.second < m_filterThreshold * getNoiseRMSPerCell(p.first)); })) != aCells.end()) { aCells.erase(it++); } @@ -142,8 +142,8 @@ StatusCode NoiseCaloCellsFromFileTool::initNoiseFromFile() { for (unsigned i = 0; i < m_numRadialLayers; i++) { elecNoiseLayerHistoName = m_elecNoiseHistoName + std::to_string(i + 1); debug() << "Getting histogram with a name " << elecNoiseLayerHistoName << endmsg; - m_histoElecNoiseConst.push_back(*dynamic_cast(noiseFile->Get(elecNoiseLayerHistoName.c_str()))); - if (m_histoElecNoiseConst.at(i).GetNbinsX() < 1) { + m_histoElecNoiseRMS.push_back(*dynamic_cast(noiseFile->Get(elecNoiseLayerHistoName.c_str()))); + if (m_histoElecNoiseRMS.at(i).GetNbinsX() < 1) { error() << "Histogram " << elecNoiseLayerHistoName << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; return StatusCode::FAILURE; @@ -151,8 +151,8 @@ StatusCode NoiseCaloCellsFromFileTool::initNoiseFromFile() { if (m_addPileup) { pileupLayerHistoName = m_pileupHistoName + std::to_string(i + 1); debug() << "Getting histogram with a name " << pileupLayerHistoName << endmsg; - m_histoPileupConst.push_back(*dynamic_cast(noiseFile->Get(pileupLayerHistoName.c_str()))); - if (m_histoPileupConst.at(i).GetNbinsX() < 1) { + m_histoPileupNoiseRMS.push_back(*dynamic_cast(noiseFile->Get(pileupLayerHistoName.c_str()))); + if (m_histoPileupNoiseRMS.at(i).GetNbinsX() < 1) { error() << "Histogram " << pileupLayerHistoName << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; return StatusCode::FAILURE; @@ -163,12 +163,12 @@ StatusCode NoiseCaloCellsFromFileTool::initNoiseFromFile() { noiseFile->Close(); // Check if we have same number of histograms (all layers) for pileup and electronics noise - if (m_histoElecNoiseConst.size() == 0) { + if (m_histoElecNoiseRMS.size() == 0) { error() << "No histograms with noise found!!!!" << endmsg; return StatusCode::FAILURE; } if (m_addPileup) { - if (m_histoElecNoiseConst.size() != m_histoPileupConst.size()) { + if (m_histoElecNoiseRMS.size() != m_histoPileupNoiseRMS.size()) { error() << "Missing histograms! Different number of histograms for electronics noise and pileup!!!!" << endmsg; return StatusCode::FAILURE; } @@ -177,14 +177,14 @@ StatusCode NoiseCaloCellsFromFileTool::initNoiseFromFile() { return StatusCode::SUCCESS; } -double NoiseCaloCellsFromFileTool::getNoiseConstantPerCell(int64_t aCellId) { +double NoiseCaloCellsFromFileTool::getNoiseRMSPerCell(uint64_t aCellId) { const dd4hep::DDSegmentation::FCCSWGridPhiEta_k4geo* segmentation = m_segmentationPhiEta; if (segmentation == nullptr) { segmentation = dynamic_cast(&m_segmentationMulti->subsegmentation(aCellId)); } - double elecNoise = 0.; - double pileupNoise = 0.; + double elecNoiseRMS = 0.; + double pileupNoiseRMS = 0.; double cellEta; if (m_useSeg) @@ -196,35 +196,35 @@ double NoiseCaloCellsFromFileTool::getNoiseConstantPerCell(int64_t aCellId) { // 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 ibin = m_histoElecNoiseConst.at(index).FindFixBin(fabs(cellEta)); + if (m_histoElecNoiseRMS.size() != 0) { + int ibin = m_histoElecNoiseRMS.at(index).FindFixBin(fabs(cellEta)); // Check that there are not more layers than the constants are provided for - if (cellLayer < m_histoElecNoiseConst.size()) { - elecNoise = m_histoElecNoiseConst.at(cellLayer).GetBinContent(ibin); + if (cellLayer < m_histoElecNoiseRMS.size()) { + elecNoiseRMS = m_histoElecNoiseRMS.at(cellLayer).GetBinContent(ibin); if (m_addPileup) { - pileupNoise = m_histoPileupConst.at(cellLayer).GetBinContent(ibin); + pileupNoiseRMS = m_histoPileupNoiseRMS.at(cellLayer).GetBinContent(ibin); } } else { - debug() + error() << "More radial layers than we have noise for!!!! Using the last layer for all histograms outside the range." << endmsg; } } else { - debug() << "No histograms with noise constants!!!!! " << endmsg; + error() << "No histograms with noise constants!!!!! " << endmsg; } // Total noise: electronics noise + pileup - double totalNoise = 0; + double totalNoiseRMS = 0; if (m_addPileup) { - totalNoise = sqrt(elecNoise*elecNoise + pileupNoise*pileupNoise) * m_scaleFactor; + totalNoiseRMS = sqrt(elecNoiseRMS*elecNoiseRMS + pileupNoiseRMS*pileupNoiseRMS) * m_scaleFactor; } else { // avoid useless math operations if no pileup - totalNoise = elecNoise * m_scaleFactor; + totalNoiseRMS = elecNoiseRMS * m_scaleFactor; } - if (totalNoise < 1e-6) { - debug() << "Zero noise: cell eta " << cellEta << " layer " << cellLayer << " noise " << totalNoise << endmsg; + if (totalNoiseRMS < 1e-6) { + warning() << "Zero noise RMS: cell eta " << cellEta << " layer " << cellLayer << " noise " << totalNoiseRMS << endmsg; } - return totalNoise; + return totalNoiseRMS; } diff --git a/RecCalorimeter/src/components/NoiseCaloCellsFromFileTool.h b/RecCalorimeter/src/components/NoiseCaloCellsFromFileTool.h index 237b3d1e..7296416e 100644 --- a/RecCalorimeter/src/components/NoiseCaloCellsFromFileTool.h +++ b/RecCalorimeter/src/components/NoiseCaloCellsFromFileTool.h @@ -49,8 +49,8 @@ class NoiseCaloCellsFromFileTool : public GaudiTool, virtual public INoiseCaloCe /// Open file and read noise histograms in the memory StatusCode initNoiseFromFile(); - /// Find the appropriate noise constant from the histogram - double getNoiseConstantPerCell(int64_t aCellID); + /// Find the appropriate noise RMS from the histogram + double getNoiseRMSPerCell(uint64_t aCellID); private: /// Handle for tool to get cell positions @@ -84,10 +84,10 @@ class NoiseCaloCellsFromFileTool : public GaudiTool, virtual public INoiseCaloCe this, "useAbsInFilter", false, "Cell filtering condition becomes: drop cell if abs(Ecell) < filterThreshold*m_cellNoise"}; /// Number of radial layers Gaudi::Property m_numRadialLayers{this, "numRadialLayers", 3, "Number of radial layers"}; - /// Histograms with pileup constants (index in array - radial layer) - std::vector m_histoPileupConst; - /// Histograms with electronics noise constants (index in array - radial layer) - std::vector m_histoElecNoiseConst; + /// Histograms with pileup RMS (index in array - radial layer) + std::vector m_histoPileupNoiseRMS; + /// Histograms with electronics noise RMS (index in array - radial layer) + std::vector m_histoElecNoiseRMS; /// Random Number Service IRndmGenSvc* m_randSvc; diff --git a/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp b/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp index a5f50379..6e5df2cd 100644 --- a/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp +++ b/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp @@ -69,8 +69,8 @@ StatusCode ReadNoiseFromFileTool::initNoiseFromFile() { error() << "File path: " << m_noiseFileName.value() << endmsg; return StatusCode::FAILURE; } - std::unique_ptr file(TFile::Open(m_noiseFileName.value().c_str(), "READ")); - if (file->IsZombie()) { + std::unique_ptr noiseFile(TFile::Open(m_noiseFileName.value().c_str(), "READ")); + if (noiseFile->IsZombie()) { error() << "Unable to open the file with the noise values!" << endmsg; error() << "File path: " << m_noiseFileName.value() << endmsg; return StatusCode::FAILURE; @@ -85,8 +85,8 @@ StatusCode ReadNoiseFromFileTool::initNoiseFromFile() { for (unsigned i = 0; i < m_numRadialLayers; i++) { elecNoiseLayerHistoName = m_elecNoiseHistoName + std::to_string(i + 1); debug() << "Getting histogram with a name " << elecNoiseLayerHistoName << endmsg; - m_histoElecNoiseConst.push_back(*dynamic_cast(file->Get(elecNoiseLayerHistoName.c_str()))); - if (m_histoElecNoiseConst.at(i).GetNbinsX() < 1) { + m_histoElecNoiseRMS.push_back(*dynamic_cast(noiseFile->Get(elecNoiseLayerHistoName.c_str()))); + if (m_histoElecNoiseRMS.at(i).GetNbinsX() < 1) { error() << "Histogram " << elecNoiseLayerHistoName << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; return StatusCode::FAILURE; @@ -94,7 +94,7 @@ StatusCode ReadNoiseFromFileTool::initNoiseFromFile() { if (m_setNoiseOffset){ elecNoiseOffsetLayerHistoName = m_elecNoiseOffsetHistoName + std::to_string(i + 1); debug() << "Getting histogram with a name " << elecNoiseOffsetLayerHistoName << endmsg; - m_histoElecNoiseOffset.push_back(*dynamic_cast(file->Get(elecNoiseOffsetLayerHistoName.c_str()))); + m_histoElecNoiseOffset.push_back(*dynamic_cast(noiseFile->Get(elecNoiseOffsetLayerHistoName.c_str()))); if (m_histoElecNoiseOffset.at(i).GetNbinsX() < 1) { error() << "Histogram " << elecNoiseOffsetLayerHistoName << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; @@ -104,8 +104,8 @@ StatusCode ReadNoiseFromFileTool::initNoiseFromFile() { if (m_addPileup) { pileupLayerHistoName = m_pileupHistoName + std::to_string(i + 1); debug() << "Getting histogram with a name " << pileupLayerHistoName << endmsg; - m_histoPileupConst.push_back(*dynamic_cast(file->Get(pileupLayerHistoName.c_str()))); - if (m_histoPileupConst.at(i).GetNbinsX() < 1) { + m_histoPileupNoiseRMS.push_back(*dynamic_cast(noiseFile->Get(pileupLayerHistoName.c_str()))); + if (m_histoPileupNoiseRMS.at(i).GetNbinsX() < 1) { error() << "Histogram " << pileupLayerHistoName << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; return StatusCode::FAILURE; @@ -113,7 +113,7 @@ StatusCode ReadNoiseFromFileTool::initNoiseFromFile() { if (m_setNoiseOffset == true){ pileupOffsetLayerHistoName = m_pileupOffsetHistoName + std::to_string(i + 1); debug() << "Getting histogram with a name " << pileupOffsetLayerHistoName << endmsg; - m_histoPileupOffset.push_back(*dynamic_cast(file->Get(pileupOffsetLayerHistoName.c_str()))); + m_histoPileupOffset.push_back(*dynamic_cast(noiseFile->Get(pileupOffsetLayerHistoName.c_str()))); if (m_histoElecNoiseOffset.at(i).GetNbinsX() < 1) { error() << "Histogram " << pileupOffsetLayerHistoName << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; @@ -123,12 +123,12 @@ StatusCode ReadNoiseFromFileTool::initNoiseFromFile() { } } // Check if we have same number of histograms (all layers) for pileup and electronics noise - if (m_histoElecNoiseConst.size() == 0 ) { + if (m_histoElecNoiseRMS.size() == 0) { error() << "No histograms with noise found!!!!" << endmsg; return StatusCode::FAILURE; } if (m_addPileup) { - if (m_histoElecNoiseConst.size() != m_histoPileupConst.size()) { + if (m_histoElecNoiseRMS.size() != m_histoPileupNoiseRMS.size()) { error() << "Missing histograms! Different number of histograms for electronics noise and pileup!!!!" << endmsg; return StatusCode::FAILURE; } @@ -137,10 +137,10 @@ StatusCode ReadNoiseFromFileTool::initNoiseFromFile() { return StatusCode::SUCCESS; } -double ReadNoiseFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) { +double ReadNoiseFromFileTool::getNoiseRMSPerCell(uint64_t aCellId) { - double elecNoise = 0.; - double pileupNoise = 0.; + double elecNoiseRMS = 0.; + double pileupNoiseRMS = 0.; // Get cell coordinates: eta and radial layer dd4hep::DDSegmentation::CellID cID = aCellId; @@ -151,11 +151,11 @@ double ReadNoiseFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) { // 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(); + if (m_histoElecNoiseRMS.size() != 0) { + int Nbins = m_histoElecNoiseRMS.at(index).GetNbinsX(); double deltaEtaBin = - (m_histoElecNoiseConst.at(index).GetBinLowEdge(Nbins) + m_histoElecNoiseConst.at(index).GetBinWidth(Nbins) - - m_histoElecNoiseConst.at(index).GetBinLowEdge(1)) / + (m_histoElecNoiseRMS.at(index).GetBinLowEdge(Nbins) + m_histoElecNoiseRMS.at(index).GetBinWidth(Nbins) - + m_histoElecNoiseRMS.at(index).GetBinLowEdge(1)) / Nbins; // find the eta bin for the cell int ibin = floor(fabs(cellEta) / deltaEtaBin) + 1; @@ -165,10 +165,10 @@ double ReadNoiseFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) { ibin = Nbins; } // Check that there are not more layers than the constants are provided for - if (cellLayer < m_histoElecNoiseConst.size()) { - elecNoise = m_histoElecNoiseConst.at(cellLayer).GetBinContent(ibin); + if (cellLayer < m_histoElecNoiseRMS.size()) { + elecNoiseRMS = m_histoElecNoiseRMS.at(cellLayer).GetBinContent(ibin); if (m_addPileup) { - pileupNoise = m_histoPileupConst.at(cellLayer).GetBinContent(ibin); + pileupNoiseRMS = m_histoPileupNoiseRMS.at(cellLayer).GetBinContent(ibin); } } else { error() @@ -180,22 +180,21 @@ double ReadNoiseFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) { } // Total noise: electronics noise + pileup - double totalNoise = sqrt(pow(elecNoise, 2) + pow(pileupNoise, 2)) * m_scaleFactor; + double totalNoiseRMS = sqrt(elecNoiseRMS*elecNoiseRMS + pileupNoiseRMS*pileupNoiseRMS) * m_scaleFactor; - if (totalNoise < 1e-6) { - warning() << "Zero noise: cell eta " << cellEta << " layer " << cellLayer << " noise " << totalNoise << endmsg; + if (totalNoiseRMS < 1e-6) { + warning() << "Zero noise: cell eta " << cellEta << " layer " << cellLayer << " noise " << totalNoiseRMS << endmsg; } - return totalNoise; + return totalNoiseRMS; } double ReadNoiseFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) { - if (!m_setNoiseOffset) - return 0.; - else { - double elecNoise = 0.; - double pileupNoise = 0.; + if (!m_setNoiseOffset) return 0.; + + double elecNoiseOffset = 0.; + double pileupNoiseOffset = 0.; // Get cell coordinates: eta and radial layer dd4hep::DDSegmentation::CellID cID = aCellId; @@ -220,9 +219,9 @@ double ReadNoiseFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) { } // Check that there are not more layers than the constants are provided for if (cellLayer < m_histoElecNoiseOffset.size()) { - elecNoise = m_histoElecNoiseOffset.at(cellLayer).GetBinContent(ibin); + elecNoiseOffset = m_histoElecNoiseOffset.at(cellLayer).GetBinContent(ibin); if (m_addPileup) { - pileupNoise = m_histoPileupOffset.at(cellLayer).GetBinContent(ibin); + pileupNoiseOffset = m_histoPileupOffset.at(cellLayer).GetBinContent(ibin); } } else { error() @@ -234,12 +233,8 @@ double ReadNoiseFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) { } // Total noise: electronics noise + pileup - double totalNoise = sqrt(pow(elecNoise, 2) + pow(pileupNoise, 2)) * m_scaleFactor; - - if (totalNoise < 1e-6) { - warning() << "Zero noise: cell eta " << cellEta << " layer " << cellLayer << " noise " << totalNoise << endmsg; - } + double totalNoiseOffset = sqrt(elecNoiseOffset*elecNoiseOffset + pileupNoiseOffset*pileupNoiseOffset) * m_scaleFactor; // shouldnt the offset be summed linearly? + // No warning is printed if offset is zero because that is the usual scenario - return totalNoise; - } + return totalNoiseOffset; } diff --git a/RecCalorimeter/src/components/ReadNoiseFromFileTool.h b/RecCalorimeter/src/components/ReadNoiseFromFileTool.h index 9fa5e1cc..33f676a0 100644 --- a/RecCalorimeter/src/components/ReadNoiseFromFileTool.h +++ b/RecCalorimeter/src/components/ReadNoiseFromFileTool.h @@ -41,7 +41,7 @@ class ReadNoiseFromFileTool : public GaudiTool, virtual public INoiseConstTool { /// Open file and read noise histograms in the memory StatusCode initNoiseFromFile(); /// Find the appropriate noise constant from the histogram - double getNoiseConstantPerCell(uint64_t aCellID); + double getNoiseRMSPerCell(uint64_t aCellID); double getNoiseOffsetPerCell(uint64_t aCellID); private: @@ -75,10 +75,10 @@ class ReadNoiseFromFileTool : public GaudiTool, virtual public INoiseConstTool { /// Factor to apply to the noise values to get them in GeV if e.g. they were produced in MeV Gaudi::Property m_scaleFactor{this, "scaleFactor", 1, "Factor to apply to the noise values"}; - /// Histograms with pileup constants (index in array - radial layer) - std::vector m_histoPileupConst; - /// Histograms with electronics noise constants (index in array - radial layer) - std::vector m_histoElecNoiseConst; + /// Histograms with pileup RMS (index in array - radial layer) + std::vector m_histoPileupNoiseRMS; + /// Histograms with electronics noise RMS (index in array - radial layer) + std::vector m_histoElecNoiseRMS; /// Histograms with pileup offset (index in array - radial layer) std::vector m_histoPileupOffset; diff --git a/RecCalorimeter/src/components/TopoCaloNoisyCells.cpp b/RecCalorimeter/src/components/TopoCaloNoisyCells.cpp index 9d7af4d8..75709431 100644 --- a/RecCalorimeter/src/components/TopoCaloNoisyCells.cpp +++ b/RecCalorimeter/src/components/TopoCaloNoisyCells.cpp @@ -44,7 +44,7 @@ StatusCode TopoCaloNoisyCells::initialize() { double readNoisyCells; double readNoisyCellsOffset; tree->SetBranchAddress("cellId", &readCellId); - tree->SetBranchAddress("noiseLevel", &readNoisyCells); + tree->SetBranchAddress("noiseLevel", &readNoisyCells); // would be better to call branch noiseRMS rather than noiseLevel tree->SetBranchAddress("noiseOffset", &readNoisyCellsOffset); for (uint i = 0; i < tree->GetEntries(); i++) { tree->GetEntry(i); diff --git a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp index 7ff80423..e079fc76 100644 --- a/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp +++ b/RecFCCeeCalorimeter/src/components/CreateFCCeeCaloNoiseLevelMap.cpp @@ -175,24 +175,24 @@ StatusCode CreateFCCeeCaloNoiseLevelMap::initialize() decoder->set(cellId, "phi", iphi); decoder->set(cellId, "theta", itheta + numCells[2]); // start from the minimum existing theta cell in this layer uint64_t id = cellId; - double noise = 0.; + double noiseRMS = 0.; double noiseOffset = 0.; if (m_fieldValuesSegmented[iSys] == m_ecalBarrelSysId) { - noise = m_ecalBarrelNoiseTool->getNoiseConstantPerCell(id); + noiseRMS = m_ecalBarrelNoiseTool->getNoiseRMSPerCell(id); noiseOffset = m_ecalBarrelNoiseTool->getNoiseOffsetPerCell(id); } else if (m_fieldValuesSegmented[iSys] == m_hcalBarrelSysId) { - noise = m_hcalBarrelNoiseTool->getNoiseConstantPerCell(id); + noiseRMS = m_hcalBarrelNoiseTool->getNoiseRMSPerCell(id); noiseOffset = m_hcalBarrelNoiseTool->getNoiseOffsetPerCell(id); } else { warning() << "Unexpected system value for phi-theta readout " << m_fieldValuesSegmented[iSys] - << ", setting noise to 0.0" << endmsg; + << ", setting noise RMS and offset to 0.0" << endmsg; } - map.insert(std::pair>(id, std::make_pair(noise, noiseOffset))); + map.insert(std::pair>(id, std::make_pair(noiseRMS, noiseOffset))); } } } @@ -221,19 +221,19 @@ StatusCode CreateFCCeeCaloNoiseLevelMap::initialize() decoder->set(cellId, "module", imodule * moduleThetaSegmentation->mergedModules(ilayer)); decoder->set(cellId, "theta", numCells[2] + itheta * moduleThetaSegmentation->mergedThetaCells(ilayer)); // start from the minimum existing theta cell in this layer uint64_t id = cellId; - double noise = 0.; + double noiseRMS = 0.; double noiseOffset = 0.; if (m_fieldValuesSegmented[iSys] == m_ecalBarrelSysId) { - noise = m_ecalBarrelNoiseTool->getNoiseConstantPerCell(id); + noiseRMS = m_ecalBarrelNoiseTool->getNoiseRMSPerCell(id); noiseOffset = m_ecalBarrelNoiseTool->getNoiseOffsetPerCell(id); } else { warning() << "Unexpected system value for module-theta readout " << m_fieldValuesSegmented[iSys] - << ", setting noise to 0.0" << endmsg; + << ", setting noise RMS and offset to 0.0" << endmsg; } - map.insert(std::pair>(id, std::make_pair(noise, noiseOffset))); + map.insert(std::pair>(id, std::make_pair(noiseRMS, noiseOffset))); } } } @@ -252,15 +252,15 @@ StatusCode CreateFCCeeCaloNoiseLevelMap::initialize() outFile->cd(); TTree tree("noisyCells", "Tree with map of noise per cell"); uint64_t saveCellId; - double saveNoiseLevel; + double saveNoiseRMS; double saveNoiseOffset; tree.Branch("cellId", &saveCellId, "cellId/l"); - tree.Branch("noiseLevel", &saveNoiseLevel); + tree.Branch("noiseLevel", &saveNoiseRMS); // would be better to call it noiseRMS tree.Branch("noiseOffset", &saveNoiseOffset); for (const auto &item : map) { saveCellId = item.first; - saveNoiseLevel = item.second.first; + saveNoiseRMS = item.second.first; saveNoiseOffset = item.second.second; tree.Fill(); } diff --git a/RecFCCeeCalorimeter/src/components/NoiseCaloCellsVsThetaFromFileTool.cpp b/RecFCCeeCalorimeter/src/components/NoiseCaloCellsVsThetaFromFileTool.cpp new file mode 100644 index 00000000..7a4039a6 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/NoiseCaloCellsVsThetaFromFileTool.cpp @@ -0,0 +1,276 @@ +#include "NoiseCaloCellsVsThetaFromFileTool.h" + +// k4geo +#include "detectorCommon/DetUtils_k4geo.h" + +// k4FWCore +#include "k4Interface/IGeoSvc.h" + +// DD4hep +#include "DD4hep/Detector.h" + +// ROOT +#include "TSystem.h" +#include "TFile.h" +#include "TH1F.h" +#include "TMath.h" + +DECLARE_COMPONENT(NoiseCaloCellsVsThetaFromFileTool) + +NoiseCaloCellsVsThetaFromFileTool::NoiseCaloCellsVsThetaFromFileTool(const std::string& type, const std::string& name, + const IInterface* parent) + : GaudiTool(type, name, parent), m_geoSvc("GeoSvc", name) { + declareInterface(this); + declareInterface(this); + declareProperty("cellPositionsTool", m_cellPositionsTool, "Handle for tool to retrieve cell positions"); +} + +StatusCode NoiseCaloCellsVsThetaFromFileTool::initialize() { + // Get 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; + } + + // Retrieve positioning tool + if (!m_cellPositionsTool.retrieve()) { + error() << "Unable to retrieve cell positions tool." << endmsg; + return StatusCode::FAILURE; + } + + // Initialize random service + if (service("RndmGenSvc", m_randSvc).isFailure()) { + error() << "Couldn't get RndmGenSvc!!!!" << endmsg; + return StatusCode::FAILURE; + } + + // Initialize gaussian random number generator + if (m_gauss.initialize(m_randSvc, Rndm::Gauss(0., 1.)).isFailure()) { + error() << "Couldn't initialize RndmGenSvc!!!!" << endmsg; + return StatusCode::FAILURE; + } + + // open and check file, read the histograms with noise constants + if (initNoiseFromFile().isFailure()) { + error() << "Couldn't open file with noise constants!!!" << endmsg; + return StatusCode::FAILURE; + } + + // Get decoder and index of layer field + // put in try... block + m_decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder(); + m_index_activeField = m_decoder->index(m_activeFieldName); + + debug() << "Filter noise threshold: " << m_filterThreshold << "*sigma" << endmsg; + + StatusCode sc = GaudiTool::initialize(); + if (sc.isFailure()) return sc; + + return sc; +} + +void NoiseCaloCellsVsThetaFromFileTool::addRandomCellNoise(std::unordered_map& aCells) { + std::for_each(aCells.begin(), aCells.end(), [this](std::pair& p) { + p.second += getNoiseOffsetPerCell(p.first); + p.second += (getNoiseRMSPerCell(p.first) * m_gauss.shoot()); + }); +} + +void NoiseCaloCellsVsThetaFromFileTool::filterCellNoise(std::unordered_map& aCells) { + // Erase a cell if it has energy below a threshold from the vector + auto it = aCells.begin(); + while ((it = std::find_if(it, aCells.end(), [this](std::pair& p) { + return m_useAbsInFilter ? bool(std::abs(p.second-getNoiseOffsetPerCell(p.first)) < m_filterThreshold * getNoiseRMSPerCell(p.first)) : bool(p.second < getNoiseOffsetPerCell(p.first) + m_filterThreshold * getNoiseRMSPerCell(p.first)); + })) != aCells.end()) { + aCells.erase(it++); + } +} + +StatusCode NoiseCaloCellsVsThetaFromFileTool::finalize() { + StatusCode sc = GaudiTool::finalize(); + return sc; +} + +StatusCode NoiseCaloCellsVsThetaFromFileTool::initNoiseFromFile() { + // Check if file exists + if (m_noiseFileName.empty()) { + error() << "Name of the file with the noise values not provided!" << endmsg; + return StatusCode::FAILURE; + } + if (gSystem->AccessPathName(m_noiseFileName.value().c_str())) { + error() << "Provided file with the noise values not found!" << endmsg; + error() << "File path: " << m_noiseFileName.value() << endmsg; + return StatusCode::FAILURE; + } + std::unique_ptr noiseFile(TFile::Open(m_noiseFileName.value().c_str(), "READ")); + if (noiseFile->IsZombie()) { + error() << "Unable to open the file with the noise values!" << endmsg; + error() << "File path: " << m_noiseFileName.value() << endmsg; + return StatusCode::FAILURE; + } else { + info() << "Using the following file with noise values: " + << m_noiseFileName.value() << endmsg; + } + + std::string elecNoiseRMSLayerHistoName, pileupNoiseRMSLayerHistoName; + std::string elecNoiseOffsetLayerHistoName, pileupNoiseOffsetLayerHistoName; + // Read the histograms with electronics noise and pileup from the file + for (unsigned i = 0; i < m_numRadialLayers; i++) { + elecNoiseRMSLayerHistoName = m_elecNoiseRMSHistoName + std::to_string(i + 1); + debug() << "Getting histogram with a name " << elecNoiseRMSLayerHistoName << endmsg; + m_histoElecNoiseRMS.push_back(*dynamic_cast(noiseFile->Get(elecNoiseRMSLayerHistoName.c_str()))); + if (m_histoElecNoiseRMS.at(i).GetNbinsX() < 1) { + error() << "Histogram " << elecNoiseRMSLayerHistoName + << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; + return StatusCode::FAILURE; + } + if (m_setNoiseOffset){ + elecNoiseOffsetLayerHistoName = m_elecNoiseOffsetHistoName + std::to_string(i + 1); + debug() << "Getting histogram with a name " << elecNoiseOffsetLayerHistoName << endmsg; + m_histoElecNoiseOffset.push_back(*dynamic_cast(noiseFile->Get(elecNoiseOffsetLayerHistoName.c_str()))); + if (m_histoElecNoiseOffset.at(i).GetNbinsX() < 1) { + error() << "Histogram " << elecNoiseOffsetLayerHistoName + << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; + return StatusCode::FAILURE; + } + } + if (m_addPileup) { + pileupNoiseRMSLayerHistoName = m_pileupNoiseRMSHistoName + std::to_string(i + 1); + debug() << "Getting histogram with a name " << pileupNoiseRMSLayerHistoName << endmsg; + m_histoPileupNoiseRMS.push_back(*dynamic_cast(noiseFile->Get(pileupNoiseRMSLayerHistoName.c_str()))); + if (m_histoPileupNoiseRMS.at(i).GetNbinsX() < 1) { + error() << "Histogram " << pileupNoiseRMSLayerHistoName + << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; + return StatusCode::FAILURE; + } + if (m_setNoiseOffset){ + pileupNoiseOffsetLayerHistoName = m_pileupNoiseOffsetHistoName + std::to_string(i + 1); + debug() << "Getting histogram with a name " << pileupNoiseOffsetLayerHistoName << endmsg; + m_histoPileupNoiseOffset.push_back(*dynamic_cast(noiseFile->Get(pileupNoiseOffsetLayerHistoName.c_str()))); + if (m_histoElecNoiseOffset.at(i).GetNbinsX() < 1) { + error() << "Histogram " << pileupNoiseOffsetLayerHistoName + << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; + return StatusCode::FAILURE; + } + } + } + } + noiseFile->Close(); + + // Check if we have same number of histograms (all layers) for pileup and electronics noise + if (m_histoElecNoiseRMS.size() == 0) { + error() << "No histograms with noise found!!!!" << endmsg; + return StatusCode::FAILURE; + } + if (m_addPileup) { + if (m_histoElecNoiseRMS.size() != m_histoPileupNoiseRMS.size()) { + error() << "Missing histograms! Different number of histograms for electronics noise RMS and pileup noise RMS!!!!" << endmsg; + return StatusCode::FAILURE; + } + } + if (m_setNoiseOffset) { + if (m_histoElecNoiseOffset.size() != m_histoElecNoiseRMS.size()) { + error() << "Missing histograms! Different number of histograms for electronics noise RMS and offset!!!!" << endmsg; + return StatusCode::FAILURE; + } + if (m_addPileup) { + if (m_histoPileupNoiseOffset.size() != m_histoElecNoiseRMS.size()) { + error() << "Missing histograms! Different number of histograms for electronics noise RMS and pileup noise offset!!!!" << endmsg; + return StatusCode::FAILURE; + } + } + } + + return StatusCode::SUCCESS; +} + +double NoiseCaloCellsVsThetaFromFileTool::getNoiseRMSPerCell(uint64_t aCellId) { + + double elecNoiseRMS = 0.; + double pileupNoiseRMS = 0.; + + double cellTheta = m_cellPositionsTool->xyzPosition(aCellId).Theta(); + unsigned cellLayer = m_decoder->get(aCellId, m_index_activeField); + + // All histograms have same binning, all bins with same size + // Using the histogram in the first layer to get the bin size + if (m_histoElecNoiseRMS.size() != 0) { + unsigned index = 0; + int ibin = m_histoElecNoiseRMS.at(index).FindBin(cellTheta); + // Check that there are not more layers than the constants are provided for + if (cellLayer < m_histoElecNoiseRMS.size()) { + elecNoiseRMS = m_histoElecNoiseRMS.at(cellLayer).GetBinContent(ibin); + if (m_addPileup) { + pileupNoiseRMS = m_histoPileupNoiseRMS.at(cellLayer).GetBinContent(ibin); + } + } else { + error() + << "More radial layers than we have noise for!!!! Using the last layer for all histograms outside the range." + << endmsg; + } + } else { + error() << "No histograms with noise constants!!!!! " << endmsg; + } + + // Total noise: electronics noise + pileup + double totalNoiseRMS = 0; + if (m_addPileup) { + totalNoiseRMS = sqrt(elecNoiseRMS*elecNoiseRMS + pileupNoiseRMS*pileupNoiseRMS) * m_scaleFactor; + } + else { // avoid useless math operations if no pileup + totalNoiseRMS = elecNoiseRMS * m_scaleFactor; + } + + if (totalNoiseRMS < 1e-6) { + warning() << "Zero noise RMS: cell theta " << cellTheta << " layer " << cellLayer << " noise RMS " << totalNoiseRMS << endmsg; + } + + return totalNoiseRMS; +} + + +double NoiseCaloCellsVsThetaFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) { + + if (!m_setNoiseOffset) return 0.; + + double elecNoiseOffset = 0.; + double pileupNoiseOffset = 0.; + + // Get cell coordinates: theta and radial layer + double cellTheta = m_cellPositionsTool->xyzPosition(aCellId).Theta(); + unsigned cellLayer = m_decoder->get(aCellId, m_activeFieldName); + + // All histograms have same binning, all bins with same size + // Using the histogram in the first layer to get the bin size + if (m_histoElecNoiseOffset.size() != 0) { + unsigned index = 0; + int Nbins = m_histoElecNoiseOffset.at(index).GetNbinsX(); + int ibin = m_histoElecNoiseOffset.at(index).FindBin(cellTheta); + if (ibin > Nbins) { + error() << "theta outside range of the histograms! Cell theta: " + << cellTheta << " Nbins in histogram: " << Nbins << endmsg; + ibin = Nbins; + } + + // Check that there are not more layers than the constants are provided for + if (cellLayer < m_histoElecNoiseOffset.size()) { + elecNoiseOffset = m_histoElecNoiseOffset.at(cellLayer).GetBinContent(ibin); + if (m_addPileup) { + pileupNoiseOffset = m_histoPileupNoiseOffset.at(cellLayer).GetBinContent(ibin); + } + } else { + error() + << "More radial layers than we have noise for!!!! Using the last layer for all histograms outside the range." + << endmsg; + } + } else { + error() << "No histograms with noise offset!!!!! " << endmsg; + } + + // Total noise offset: electronics noise + pileup (linear or quadratic sum??) + // double totalNoiseOffset = sqrt(pow(elecNoiseOffset, 2) + pow(pileupNoiseOffset, 2)) * m_scaleFactor; + double totalNoiseOffset = (elecNoiseOffset + pileupNoiseOffset) * m_scaleFactor; + + return totalNoiseOffset; +} diff --git a/RecFCCeeCalorimeter/src/components/NoiseCaloCellsVsThetaFromFileTool.h b/RecFCCeeCalorimeter/src/components/NoiseCaloCellsVsThetaFromFileTool.h new file mode 100644 index 00000000..fd0506e1 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/NoiseCaloCellsVsThetaFromFileTool.h @@ -0,0 +1,127 @@ +#ifndef RECFCCEECALORIMETER_NOISECALOCELLSVSTHETAFROMFILETOOL_H +#define RECFCCEECALORIMETER_NOISECALOCELLSVSTHETAFROMFILETOOL_H + +// from Gaudi +#include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/RndmGenerators.h" +#include "GaudiKernel/IRndmGenSvc.h" + +// k4geo +// #include "detectorSegmentations/FCCSWGridPhiEta_k4geo.h" + +// k4FWCore +#include "k4Interface/INoiseConstTool.h" +#include "k4Interface/INoiseCaloCellsTool.h" +#include "k4Interface/ICellPositionsTool.h" + +class IGeoSvc; + +// Root +class TH1F; + +/** @class NoiseCaloCellsVsThetaFromFileTool + * + * Tool for calorimeter noise + * Access noise constants from TH1F histogram (noise vs. theta) + * createRandomCellNoise: Create random CaloHits (gaussian distribution) for the vector of cells + * filterCellNoise: remove cells with energy below threshold*sigma from the vector of cells + * The tool needs a cell positioning tool to translate cellID to cell theta. + * In alternative, the tool could be rewritten to use a specific segmentation class for the cellID->theta + * translation, but it would be coupled to a specific readout. + * To avoid all these calls to the positioning tool, one could either + * - save directly the noise histograms as histos of noise vs thetaID + * - or, keep histos of noise vs theta, but change the interfaces and the tool to accept + * cells rather than cellIDs as input. One would then get theta from the cells. + * + * @author Giovanni Marchiori + * @date 2024-07 + * + */ + +class NoiseCaloCellsVsThetaFromFileTool : public GaudiTool, virtual public INoiseCaloCellsTool, virtual public INoiseConstTool { +public: + NoiseCaloCellsVsThetaFromFileTool(const std::string& type, const std::string& name, const IInterface* parent); + virtual ~NoiseCaloCellsVsThetaFromFileTool() = default; + virtual StatusCode initialize() final; + virtual StatusCode finalize() final; + + /** @brief Create random CaloHits (gaussian distribution) for the vector of cells (aCells). + * Vector of cells must contain all cells in the calorimeter with their cellIDs. + */ + virtual void addRandomCellNoise(std::unordered_map& aCells) final; + /** @brief Remove cells with energy below threshold*sigma from the vector of cells + */ + virtual void filterCellNoise(std::unordered_map& aCells) final; + + /// Open file and read noise histograms in the memory + StatusCode initNoiseFromFile(); + /// Find the appropriate noise RMS from the histogram + double getNoiseRMSPerCell(uint64_t aCellID); + double getNoiseOffsetPerCell(uint64_t aCellID); + +private: + /// Handle for tool to get cell positions + ToolHandle m_cellPositionsTool{"CellPositionsDummyTool", this}; + + /// Add pileup contribution to the electronics noise? (only if read from file) + Gaudi::Property m_addPileup{this, "addPileup", true, + "Add pileup contribution to the electronics noise? (only if read from file)"}; + /// Read noise offset from histograms in files or not (if false, offset is set to 0) + Gaudi::Property m_setNoiseOffset{this, "setNoiseOffset", false, "Set a noise offset per cell"}; + /// Name of the file with noise constants + Gaudi::Property m_noiseFileName{this, "noiseFileName", "", "Name of the file with noise constants"}; + /// Factor to apply to the noise values to get them in GeV if e.g. they were produced in MeV + Gaudi::Property m_scaleFactor{this, "scaleFactor", 1, "Factor to apply to the noise values"}; + /// Name of the detector readout (only needed to retrieve the decoder) + Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelThetaModuleMerged", + "Name of the detector readout"}; + /// Name of active layers for sampling calorimeter + Gaudi::Property m_activeFieldName{this, "activeFieldName", "active_layer", + "Name of active layers for sampling calorimeter"}; + /// Name of pileup noise RMS histogram + Gaudi::Property m_pileupNoiseRMSHistoName{this, "pileupNoiseRMSHistoName", "h_pileup_layer", + "Name of pileup noise RMS histogram"}; + /// Name of electronics noise RMS histogram + Gaudi::Property m_elecNoiseRMSHistoName{this, "elecNoiseRMSHistoName", "h_elecNoise_layer", + "Name of electronics noise RMS histogram"}; + /// Name of electronics noise offset histogram + Gaudi::Property m_elecNoiseOffsetHistoName{this, "elecNoiseOffsetHistoName", "h_mean_pileup_layer", + "Name of electronics noise offset histogram"}; + /// Name of pileup noise offset histogram + Gaudi::Property m_pileupNoiseOffsetHistoName{this, "pileupNoiseOffsetHistoName", "h_pileup_layer", + "Name of pileup noise offset histogram"}; + + /// Energy threshold (cells with Ecell < filterThreshold*m_cellNoise removed) + Gaudi::Property m_filterThreshold{this, "filterNoiseThreshold", 3, + "Energy threshold (cells with Ecell < offset + filterThreshold*m_cellNoise removed)"}; + /// Change the cell filter condition to remove only cells with abs(Ecell) < offset + filterThreshold*m_cellNoise removed) + /// This avoids to keep only 'one side' of the noise fluctuations and prevents biasing cluster energy towards higher energies + Gaudi::Property m_useAbsInFilter{this, "useAbsInFilter", false, + "If set, cell filtering condition becomes: drop cell if abs(Ecell-offset) < filterThreshold*m_cellNoise"}; + /// Number of radial layers + Gaudi::Property m_numRadialLayers{this, "numRadialLayers", 3, + "Number of radial layers"}; + + /// Histograms with pileup noise RMS (index in array - radial layer) + std::vector m_histoPileupNoiseRMS; + /// Histograms with electronics noise RMS (index in array - radial layer) + std::vector m_histoElecNoiseRMS; + /// Histograms with pileup offset (index in array - radial layer) + std::vector m_histoPileupNoiseOffset; + /// Histograms with electronics noise offset (index in array - radial layer) + std::vector m_histoElecNoiseOffset; + + /// Random Number Service + IRndmGenSvc* m_randSvc; + /// Gaussian random number generator used for the generation of random noise hits + Rndm::Numbers m_gauss; + + /// Pointer to the geometry service + ServiceHandle m_geoSvc; + + /// Decoder for ECal layers + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + int m_index_activeField; +}; + +#endif /* RECFCCEECALORIMETER_NOISECALOCELLSVSTHETAFROMFILETOOL_H */ diff --git a/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.cpp b/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.cpp deleted file mode 100644 index d461cb83..00000000 --- a/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.cpp +++ /dev/null @@ -1,270 +0,0 @@ -#include "ReadNoiseVsThetaFromFileTool.h" - -// k4geo -#include "detectorCommon/DetUtils_k4geo.h" - -// k4FWCore -#include "k4Interface/IGeoSvc.h" - - -// DD4hep -#include "DD4hep/Detector.h" -#include "DD4hep/Readout.h" -#include "DDSegmentation/Segmentation.h" - -// ROOT -#include "TSystem.h" -#include "TFile.h" -#include "TH1F.h" - -DECLARE_COMPONENT(ReadNoiseVsThetaFromFileTool) - -ReadNoiseVsThetaFromFileTool::ReadNoiseVsThetaFromFileTool(const std::string& type, const std::string& name, const IInterface* parent) - : GaudiTool(type, name, parent) { - declareInterface(this); - declareProperty("cellPositionsTool", m_cellPositionsTool, "Handle for tool to retrieve cell positions"); -} - -StatusCode ReadNoiseVsThetaFromFileTool::initialize() { - // Get GeoSvc - 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; - } - - // Check if cell position tool available if m_useSeg==false; if tool not - // available, try using segmentation instead - if (!m_useSeg){ - if (!m_cellPositionsTool.retrieve()) { - error() << "Unable to retrieve cell positions tool, and useSegmentation is false." << endmsg; - return StatusCode::FAILURE; - } - } - - // Get segmentation - if (m_useSeg) { - m_segmentation = dynamic_cast( - m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation()); - if (m_segmentation == nullptr) { - error() << "There is no module-theta segmentation." << endmsg; - return StatusCode::FAILURE; - } - else - info() << "Found module-theta segmentation." << endmsg; - } - - // open and check file, read the histograms with noise constants - if (ReadNoiseVsThetaFromFileTool::initNoiseFromFile().isFailure()) { - error() << "Couldn't open file with noise constants!!!" << endmsg; - return StatusCode::FAILURE; - } - - // 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; - - return sc; -} - -StatusCode ReadNoiseVsThetaFromFileTool::finalize() { - StatusCode sc = GaudiTool::finalize(); - return sc; -} - -StatusCode ReadNoiseVsThetaFromFileTool::initNoiseFromFile() { - // check if file exists - if (m_noiseFileName.empty()) { - error() << "Name of the file with the noise values not provided!" << endmsg; - return StatusCode::FAILURE; - } - if (gSystem->AccessPathName(m_noiseFileName.value().c_str())) { - error() << "Provided file with the noise values not found!" << endmsg; - error() << "File path: " << m_noiseFileName.value() << endmsg; - return StatusCode::FAILURE; - } - std::unique_ptr inFile(TFile::Open(m_noiseFileName.value().c_str(), "READ")); - if (inFile->IsZombie()) { - error() << "Couldn't open the file with the noise values!" << endmsg; - error() << "File path: " << m_noiseFileName.value() << endmsg; - return StatusCode::FAILURE; - } else { - info() << "Opening the following file with the noise values: " - << m_noiseFileName << endmsg; - } - - std::string elecNoiseLayerHistoName, pileupLayerHistoName; - std::string elecNoiseOffsetLayerHistoName, pileupOffsetLayerHistoName; - // Read the histograms with electronics noise and pileup from the file - for (unsigned i = 0; i < m_numRadialLayers; i++) { - elecNoiseLayerHistoName = m_elecNoiseHistoName + std::to_string(i + 1); - debug() << "Getting histogram with a name " << elecNoiseLayerHistoName << endmsg; - m_histoElecNoiseConst.push_back(*dynamic_cast(inFile->Get(elecNoiseLayerHistoName.c_str()))); - if (m_histoElecNoiseConst.at(i).GetNbinsX() < 1) { - error() << "Histogram " << elecNoiseLayerHistoName - << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; - return StatusCode::FAILURE; - } - if (m_setNoiseOffset){ - elecNoiseOffsetLayerHistoName = m_elecNoiseOffsetHistoName + std::to_string(i + 1); - debug() << "Getting histogram with a name " << elecNoiseOffsetLayerHistoName << endmsg; - m_histoElecNoiseOffset.push_back(*dynamic_cast(inFile->Get(elecNoiseOffsetLayerHistoName.c_str()))); - if (m_histoElecNoiseOffset.at(i).GetNbinsX() < 1) { - error() << "Histogram " << elecNoiseOffsetLayerHistoName - << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; - return StatusCode::FAILURE; - } - } - if (m_addPileup) { - pileupLayerHistoName = m_pileupHistoName + std::to_string(i + 1); - debug() << "Getting histogram with a name " << pileupLayerHistoName << endmsg; - m_histoPileupConst.push_back(*dynamic_cast(inFile->Get(pileupLayerHistoName.c_str()))); - if (m_histoPileupConst.at(i).GetNbinsX() < 1) { - error() << "Histogram " << pileupLayerHistoName - << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg; - return StatusCode::FAILURE; - } - if (m_setNoiseOffset == true){ - pileupOffsetLayerHistoName = m_pileupOffsetHistoName + std::to_string(i + 1); - debug() << "Getting histogram with a name " << pileupOffsetLayerHistoName << endmsg; - m_histoPileupOffset.push_back(*dynamic_cast(inFile->Get(pileupOffsetLayerHistoName.c_str()))); - if (m_histoElecNoiseOffset.at(i).GetNbinsX() < 1) { - error() << "Histogram " << pileupOffsetLayerHistoName - << " 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_histoElecNoiseConst.size() == 0 ) { - error() << "No histograms with noise found!!!!" << endmsg; - return StatusCode::FAILURE; - } - if (m_addPileup) { - if (m_histoElecNoiseConst.size() != m_histoPileupConst.size()) { - error() << "Missing histograms! Different number of histograms for electronics noise and pileup!!!!" << endmsg; - return StatusCode::FAILURE; - } - } - - return StatusCode::SUCCESS; -} - -double ReadNoiseVsThetaFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) { - - double elecNoise = 0.; - double pileupNoise = 0.; - - // Get cell coordinates: eta/theta and radial layer - dd4hep::DDSegmentation::CellID cID = aCellId; - double cellTheta; - // checked that for baseline theta-module merged segmentation - // the two approaches give identical theta. - // however, code based on positioning tool is more general - // since it can be run for any segmentation class without the need - // to change the interface of this tool.. - if (m_useSeg) - cellTheta = m_segmentation->theta(aCellId); - else - cellTheta = m_cellPositionsTool->xyzPosition(aCellId).Theta(); - - 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(); - int ibin = m_histoElecNoiseConst.at(index).FindBin(cellTheta); - if (ibin > Nbins) { - error() << "theta outside range of the histograms! Cell theta: " << cellTheta << " Nbins in histogram: " << Nbins << endmsg; - ibin = Nbins; - } - // Check that there are not more layers than the constants are provided for - if (cellLayer < m_histoElecNoiseConst.size()) { - elecNoise = m_histoElecNoiseConst.at(cellLayer).GetBinContent(ibin); - if (m_addPileup) { - pileupNoise = m_histoPileupConst.at(cellLayer).GetBinContent(ibin); - } - } else { - error() - << "More radial layers than we have noise for!!!! Using the last layer for all histograms outside the range." - << endmsg; - } - } else { - error() << "No histograms with noise constants!!!!! " << endmsg; - } - - // Total noise: electronics noise + pileup - double totalNoise = sqrt(pow(elecNoise, 2) + pow(pileupNoise, 2)) * m_scaleFactor; - - if (totalNoise < 1e-6) { - warning() << "Zero noise: cell theta " << cellTheta << " layer " - << cellLayer << " noise " << totalNoise << endmsg; - } - - return totalNoise; -} - -double ReadNoiseVsThetaFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) { - - if (!m_setNoiseOffset) - return 0.; - else { - double elecNoise = 0.; - double pileupNoise = 0.; - - // Get cell coordinates: eta and radial layer - dd4hep::DDSegmentation::CellID cID = aCellId; - double cellTheta; - if (m_useSeg) - cellTheta = m_segmentation->theta(aCellId); - else - cellTheta = m_cellPositionsTool->xyzPosition(aCellId).Theta(); - 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(); - int ibin = m_histoElecNoiseOffset.at(index).FindBin(cellTheta); - if (ibin > Nbins) { - error() << "theta outside range of the histograms! Cell theta: " - << cellTheta << " Nbins in histogram: " << Nbins << endmsg; - ibin = Nbins; - } - - // Check that there are not more layers than the constants are provided for - if (cellLayer < m_histoElecNoiseOffset.size()) { - elecNoise = m_histoElecNoiseOffset.at(cellLayer).GetBinContent(ibin); - if (m_addPileup) { - pileupNoise = m_histoPileupOffset.at(cellLayer).GetBinContent(ibin); - } - } else { - error() - << "More radial layers than we have noise for!!!! Using the last layer for all histograms outside the range." - << endmsg; - } - } else { - error() << "No histograms with noise offset!!!!! " << endmsg; - } - - // Total noise: electronics noise + pileup - double totalNoise = sqrt(pow(elecNoise, 2) + pow(pileupNoise, 2)) * m_scaleFactor; - - if (totalNoise < 1e-6) { - warning() << "Zero noise: cell theta " << cellTheta << " layer " << cellLayer << " noise " << totalNoise << endmsg; - } - - return totalNoise; - } -} diff --git a/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.h b/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.h deleted file mode 100644 index 2e31669a..00000000 --- a/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.h +++ /dev/null @@ -1,98 +0,0 @@ -#ifndef RECFCCEECALORIMETER_READNOISEVSTHETAFROMFILETOOL_H -#define RECFCCEECALORIMETER_READNOISEVSTHETAFROMFILETOOL_H - -// Gaudi -#include "GaudiAlg/GaudiTool.h" -#include "GaudiKernel/ToolHandle.h" - -// k4FWCore -#include "k4Interface/INoiseConstTool.h" -#include "k4Interface/ICellPositionsTool.h" - -// k4geo -#include "detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h" - -class IGeoSvc; - -// Root -class TH1F; - -/** @class ReadNoiseVsThetaFromFileTool - * - * Tool to read the stored noise constant per cell in the calorimeters - * Access noise constants from TH1F histogram (noise vs. theta) - * Cell positioning can be done via segmentation or via a cell - * positioning tool - * - * @author Tong Li, Giovanni Marchiori - * @date 2023-09 - * - */ - -class ReadNoiseVsThetaFromFileTool : public GaudiTool, virtual public INoiseConstTool { -public: - ReadNoiseVsThetaFromFileTool(const std::string& type, const std::string& name, const IInterface* parent); - virtual ~ReadNoiseVsThetaFromFileTool() = default; - - virtual StatusCode initialize() final; - - virtual StatusCode finalize() final; - - /// Open file and read noise histograms in the memory - StatusCode initNoiseFromFile(); - /// Find the appropriate noise constant from the histogram - double getNoiseConstantPerCell(uint64_t aCellID); - double getNoiseOffsetPerCell(uint64_t aCellID); - -private: - /// Handle for tool to get cell positions (optional - not needed if segmentation class is used) - /// Provided as a way to use the tool for other theta-based segmentations (using the positioning - /// tool instead of the segmentation to compute the theta of a cell) - /// (maybe it suffices to use a GridTheta segmentation base class?) - ToolHandle m_cellPositionsTool; - /// use segmentation (theta-module only so far!) in case no cell position tool is defined - Gaudi::Property 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 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 m_setNoiseOffset{this, "setNoiseOffset", true, "Set a noise offset per cell"}; - /// Name of the file with noise constants - Gaudi::Property m_noiseFileName{this, "noiseFileName", "", "Name of the file with noise constants"}; - /// Name of the detector readout (needed to get the decoder to extract e.g. layer information) - Gaudi::Property m_readoutName{this, "readoutName", "ECalBarrelThetaModuleMerged", "Name of the detector readout"}; - /// Name of active layers for sampling calorimeter - Gaudi::Property m_activeFieldName{this, "activeFieldName", "active_layer", - "Name of active layers for sampling calorimeter"}; - /// Name of pileup histogram - Gaudi::Property m_pileupHistoName{this, "pileupHistoName", "h_pileup_layer", "Name of pileup histogram"}; - /// Name of electronics noise histogram - Gaudi::Property m_elecNoiseHistoName{this, "elecNoiseHistoName", "h_elecNoise_layer", - "Name of electronics noise histogram"}; - /// Name of electronics noise offset histogram - Gaudi::Property m_elecNoiseOffsetHistoName{this, "elecNoiseOffsetHistoName", "h_mean_pileup_layer", - "Name of electronics noise offset histogram"}; - /// Name of pileup offset histogram - Gaudi::Property m_pileupOffsetHistoName{this, "pileupOffsetHistoName", "h_pileup_layer", "Name of pileup offset histogram"}; - /// Number of radial layers - Gaudi::Property 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 m_scaleFactor{this, "scaleFactor", 1, "Factor to apply to the noise values"}; - /// Histograms with pileup constants (index in array - radial layer) - std::vector m_histoPileupConst; - /// Histograms with electronics noise constants (index in array - radial layer) - std::vector m_histoElecNoiseConst; - /// Histograms with pileup offset (index in array - radial layer) - std::vector m_histoPileupOffset; - /// Histograms with electronics noise offset (index in array - radial layer) - std::vector m_histoElecNoiseOffset; - /// Pointer to the geometry service - SmartIF m_geoSvc; - /// Theta-Module segmentation - dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo* m_segmentation; - // Decoder for ECal layers - dd4hep::DDSegmentation::BitFieldCoder* m_decoder; - -}; - -#endif /* RECFCCEECALORIMETER_READNOISEVSTHETAFROMFILETOOL_H */ diff --git a/RecFCChhCalorimeter/src/components/CreateFCChhCaloNoiseLevelMap.cpp b/RecFCChhCalorimeter/src/components/CreateFCChhCaloNoiseLevelMap.cpp index fa587faf..9552d108 100644 --- a/RecFCChhCalorimeter/src/components/CreateFCChhCaloNoiseLevelMap.cpp +++ b/RecFCChhCalorimeter/src/components/CreateFCChhCaloNoiseLevelMap.cpp @@ -90,16 +90,16 @@ StatusCode CreateFCChhCaloNoiseLevelMap::initialize() { decoder->set(cellId, "phi", iphi); decoder->set(cellId, "eta", ieta + numCells[2]); // start from the minimum existing eta cell in this layer uint64_t id = cellId; - double noise = 0.; + double noiseRMS = 0.; double noiseOffset = 0.; if (m_fieldValuesSegmented[iSys] == m_hcalBarrelSysId){ - noise = m_hcalBarrelNoiseTool->getNoiseConstantPerCell(id); + noiseRMS = m_hcalBarrelNoiseTool->getNoiseRMSPerCell(id); noiseOffset = m_hcalBarrelNoiseTool->getNoiseOffsetPerCell(id); } else if (m_fieldValuesSegmented[iSys] == m_ecalBarrelSysId){ - noise = m_ecalBarrelNoiseTool->getNoiseConstantPerCell(id); + noiseRMS = m_ecalBarrelNoiseTool->getNoiseRMSPerCell(id); noiseOffset = m_ecalBarrelNoiseTool->getNoiseOffsetPerCell(id); } - map.insert( std::pair >(id, std::make_pair(noise, noiseOffset) ) ); + map.insert( std::pair >(id, std::make_pair(noiseRMS, noiseOffset) ) ); } } } @@ -190,10 +190,10 @@ StatusCode CreateFCChhCaloNoiseLevelMap::initialize() { decoder->set(cID, m_activeFieldNamesNested[1], iphi); decoder->set(cID, m_activeFieldNamesNested[2], iz); - double noise = m_hcalBarrelNoiseTool->getNoiseConstantPerCell(cID); + double noiseRMS = m_hcalBarrelNoiseTool->getNoiseRMSPerCell(cID); double noiseOffset = m_hcalBarrelNoiseTool->getNoiseOffsetPerCell(cID); - map.insert( std::pair >(cID, std::make_pair(noise, noiseOffset) ) ); + map.insert( std::pair >(cID, std::make_pair(noiseRMS, noiseOffset) ) ); } } }