From 68792c4db73417228b50acb52c85f2217c177438 Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori Date: Fri, 29 Sep 2023 16:49:52 +0200 Subject: [PATCH] move ReadNoiseFromFileTool for theta-based readout to new tool in FCCee subdirectory --- .../src/components/ReadNoiseFromFileTool.cpp | 97 ++----- .../src/components/ReadNoiseFromFileTool.h | 19 +- .../src/components/CaloTopoClusterFCCee.h | 3 +- .../ReadNoiseVsThetaFromFileTool.cpp | 259 ++++++++++++++++++ .../components/ReadNoiseVsThetaFromFileTool.h | 99 +++++++ 5 files changed, 388 insertions(+), 89 deletions(-) create mode 100644 RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.cpp create mode 100644 RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.h diff --git a/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp b/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp index 379b2025..032b88b9 100644 --- a/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp +++ b/RecCalorimeter/src/components/ReadNoiseFromFileTool.cpp @@ -19,7 +19,6 @@ DECLARE_COMPONENT(ReadNoiseFromFileTool) ReadNoiseFromFileTool::ReadNoiseFromFileTool(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 ReadNoiseFromFileTool::initialize() { @@ -30,26 +29,10 @@ StatusCode ReadNoiseFromFileTool::initialize() { << "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()) { - info() << "Unable to retrieve cell positions tool, try eta-phi segmentation." << endmsg; - m_useSeg = true; - } - } - - // Get PhiEta segmentation - if (m_useSeg) { - m_segmentation = dynamic_cast( - m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation()); - if (m_segmentation == nullptr) { - error() << "There is no phi-eta segmentation." << endmsg; - return StatusCode::FAILURE; - } - else - info() << "Found phi-eta segmentation." << endmsg; + m_segmentation = dynamic_cast(m_geoSvc->getDetector()->readout(m_readoutName).segmentation().segmentation()); + if (m_segmentation == nullptr) { + error() << "There is no phi-eta segmentation!!!!" << endmsg; + return StatusCode::FAILURE; } // open and check file, read the histograms with noise constants @@ -60,10 +43,6 @@ StatusCode ReadNoiseFromFileTool::initialize() { // Take readout bitfield decoder from GeoSvc m_decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder(); - if (m_decoder == nullptr) { - error() << "Cannot create decore for readout " << m_readoutName << endmsg; - return StatusCode::FAILURE; - } StatusCode sc = GaudiTool::initialize(); if (sc.isFailure()) return sc; @@ -153,13 +132,10 @@ double ReadNoiseFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) { double elecNoise = 0.; double pileupNoise = 0.; - // Get cell coordinates: eta/theta and radial layer + // Get cell coordinates: eta and radial layer dd4hep::DDSegmentation::CellID cID = aCellId; - double cellEta, cellTheta; - if (m_useSeg) - cellEta = m_segmentation->eta(aCellId); - else - cellTheta = m_cellPositionsTool->xyzPosition(aCellId).Theta(); + double cellEta = m_segmentation->eta(cID); + unsigned cellLayer = m_decoder->get(cID, m_activeFieldName); // All histograms have same binning, all bins with same size @@ -167,31 +143,16 @@ double ReadNoiseFromFileTool::getNoiseConstantPerCell(uint64_t aCellId) { unsigned index = 0; if (m_histoElecNoiseConst.size() != 0) { int Nbins = m_histoElecNoiseConst.at(index).GetNbinsX(); - // GM: basically the same as FindBin ... - /* double deltaEtaBin = (m_histoElecNoiseConst.at(index).GetBinLowEdge(Nbins) + m_histoElecNoiseConst.at(index).GetBinWidth(Nbins) - m_histoElecNoiseConst.at(index).GetBinLowEdge(1)) / Nbins; // find the eta bin for the cell int ibin = floor(fabs(cellEta) / deltaEtaBin) + 1; - */ - int ibin; - if (m_useSeg) { - ibin = m_histoElecNoiseConst.at(index).FindBin(fabs(cellEta)); - if (ibin > Nbins) { - error() << "eta outside range of the histograms! Cell eta: " << cellEta << " Nbins in histogram: " << Nbins - << endmsg; - ibin = Nbins; - } - } - else { - 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; - } + if (ibin > Nbins) { + error() << "eta outside range of the histograms! Cell eta: " << cellEta << " 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()) { @@ -228,11 +189,7 @@ double ReadNoiseFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) { // Get cell coordinates: eta and radial layer dd4hep::DDSegmentation::CellID cID = aCellId; - double cellEta, cellTheta; - if (m_useSeg) - cellEta = m_segmentation->eta(aCellId); - else - cellTheta = m_cellPositionsTool->xyzPosition(aCellId).Theta(); + double cellEta = m_segmentation->eta(cID); unsigned cellLayer = m_decoder->get(cID, m_activeFieldName); // All histograms have same binning, all bins with same size @@ -240,33 +197,17 @@ double ReadNoiseFromFileTool::getNoiseOffsetPerCell(uint64_t aCellId) { unsigned index = 0; if (m_histoElecNoiseOffset.size() != 0) { int Nbins = m_histoElecNoiseOffset.at(index).GetNbinsX(); - // find the eta bin for the cell - // GM: basically the same as FindBin ... - /* - double deltaEtaBin = + double deltaEtaBin = (m_histoElecNoiseOffset.at(index).GetBinLowEdge(Nbins) + m_histoElecNoiseOffset.at(index).GetBinWidth(Nbins) - m_histoElecNoiseOffset.at(index).GetBinLowEdge(1)) / Nbins; - int ibin = floor(fabs(cellEta) / deltaEtaBin) + 1; - */ - int ibin; - if (m_useSeg) { - ibin = m_histoElecNoiseOffset.at(index).FindBin(fabs(cellEta)); - if (ibin > Nbins) { - error() << "eta outside range of the histograms! Cell eta: " << cellEta << " Nbins in histogram: " << Nbins - << endmsg; - ibin = Nbins; - } - } - else { - 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; - } + // find the eta bin for the cell + int ibin = floor(fabs(cellEta) / deltaEtaBin) + 1; + if (ibin > Nbins) { + error() << "eta outside range of the histograms! Cell eta: " << cellEta << " 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); diff --git a/RecCalorimeter/src/components/ReadNoiseFromFileTool.h b/RecCalorimeter/src/components/ReadNoiseFromFileTool.h index 189eb898..fbbc55b0 100644 --- a/RecCalorimeter/src/components/ReadNoiseFromFileTool.h +++ b/RecCalorimeter/src/components/ReadNoiseFromFileTool.h @@ -7,12 +7,11 @@ // FCCSW #include "k4FWCore/DataHandle.h" +#include "DetSegmentation/FCCSWGridPhiEta.h" #include "k4Interface/ICalorimeterTool.h" #include "k4Interface/INoiseConstTool.h" #include "k4Interface/ICellPositionsTool.h" -#include "DetSegmentation/FCCSWGridPhiEta.h" - class IGeoSvc; // Root @@ -21,7 +20,7 @@ class TH1F; /** @class ReadNoiseFromFileTool * * Tool to read the stored noise constant per cell in the calorimeters - * Access noise constants from TH1F histogram (noise vs. |eta| or theta) + * Access noise constants from TH1F histogram (noise vs. |eta|) * * @author Jana Faltova, Coralie Neubueser * @date 2018-01 @@ -44,18 +43,15 @@ class ReadNoiseFromFileTool : public GaudiTool, virtual public INoiseConstTool { double getNoiseOffsetPerCell(uint64_t aCellID); private: - /// Handle for tool to get cell positions - ToolHandle m_cellPositionsTool; - /// use segmentation (eta-phi 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) + /// Name of the detector readout Gaudi::Property m_readoutName{this, "readoutName", "ECalHitsPhiEta", "Name of the detector readout"}; /// Name of active layers for sampling calorimeter Gaudi::Property m_activeFieldName{this, "activeFieldName", "active_layer", @@ -70,23 +66,28 @@ class ReadNoiseFromFileTool : public GaudiTool, virtual public INoiseConstTool { "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; /// PhiEta segmentation dd4hep::DDSegmentation::FCCSWGridPhiEta* m_segmentation; - // Decoder for ECal layers + // Decoder dd4hep::DDSegmentation::BitFieldCoder* m_decoder; }; diff --git a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h index fa46358b..52a07d44 100644 --- a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h +++ b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h @@ -31,8 +31,7 @@ class Segmentation; } } -/** @class CaloTopoClusterFCCeeAlgorithm Reconstruction/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h - * CombinedCaloTopoCluster.h +/** @class CaloTopoClusterFCCee k4RecCalorimeter/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h * * Algorithm building the topological clusters for the energy reconstruction, following ATLAS note * ATL-LARG-PUB-2008-002. diff --git a/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.cpp b/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.cpp new file mode 100644 index 00000000..3ae2c2e7 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.cpp @@ -0,0 +1,259 @@ +#include "ReadNoiseVsThetaFromFileTool.h" + +// FCCSW +#include "DetCommon/DetUtils.h" +#include "k4Interface/IGeoSvc.h" +#include "DDSegmentation/Segmentation.h" + +// DD4hep +#include "DD4hep/Detector.h" +#include "DD4hep/Readout.h" + +// Root +#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 noise values not set" << endmsg; + return StatusCode::FAILURE; + } + std::unique_ptr file(TFile::Open(m_noiseFileName.value().c_str(), "READ")); + if (file->IsZombie()) { + error() << "Couldn't open the file with noise constants" << endmsg; + return StatusCode::FAILURE; + } else { + info() << "Opening the file with noise constants: " << 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(file->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(file->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(file->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(file->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 new file mode 100644 index 00000000..c39c4cb5 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/ReadNoiseVsThetaFromFileTool.h @@ -0,0 +1,99 @@ +#ifndef RECFCCEECALORIMETER_READNOISEVSTHETAFROMFILETOOL_H +#define RECFCCEECALORIMETER_READNOISEVSTHETAFROMFILETOOL_H + +// from Gaudi +#include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/ToolHandle.h" + +// FCCSW +//#include "k4FWCore/DataHandle.h" +//#include "k4Interface/ICalorimeterTool.h" +#include "k4Interface/INoiseConstTool.h" +#include "k4Interface/ICellPositionsTool.h" + +#include "DetSegmentation/FCCSWGridModuleThetaMerged.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* m_segmentation; + // Decoder for ECal layers + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + +}; + +#endif /* RECFCCEECALORIMETER_READNOISEVSTHETAFROMFILETOOL_H */