Consolidate tools handling noise (#98)

* make the code agnostic of detector IDs, remove unneeded dependency on specific segmentation class, improve naming of some methods

* Const -> NoiseRMS

* add noise offset, improve name of members

* Const(ant) -> (Noise)RMS

* Const(ant) -> (Noise)RMS

* Noise(Constant) -> NoiseRMS

* NoiseConst => NoiseRMS

* improve readability

* remove duplicate tool

* set noise offset off by default. Also, add setNoiseOffset flag to ConstNoiseTool and set it to false by default to avoid having to specify the list of offsets in python

* add a comment

* implement comments from review

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; iDet<m_detectors.size(); iDet++) {
+    std::string detName = "DetID_" + m_detectors[iDet];
+    try {
+      int detID = m_geoSvc->getDetector()->constant<int>(detName);
+      m_systemNoiseRMSMap.emplace(detID, m_detectorsNoiseRMS[iDet]);
+      debug() << "Set noise RMS for detector " << detName << " (ID=" << detID << ") to: "
+	      << m_detectorsNoiseRMS[iDet] << endmsg;
+      if (m_setNoiseOffset) {
+	m_systemNoiseOffsetMap.emplace(detID, m_detectorsNoiseOffset[iDet]);
+	debug() << "Set noise offset for detector " << detName << " (ID=" << detID << ") to: "
+	        << m_detectorsNoiseOffset[iDet] << endmsg;
+      }
+    } catch (const std::exception& e) {
+      error() << "Error: detector with name " << detName << " not found, exception raised: " << e.what() << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+
   StatusCode sc = GaudiTool::initialize();
   if (sc.isFailure()) return sc;
 
@@ -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;
 }

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<uint,double> m_systemNoiseConstMap;
+  std::map<uint,double> m_systemNoiseRMSMap;
   std::map<uint,double> m_systemNoiseOffsetMap;
 
-  /// Cell thresholds (1Sigma) in GeV:
+  /// List of subdetector names (they must match what is defined in DectDimension)
+ Gaudi::Property<std::vector<std::string>> m_detectors{this, "detectors", {"ECAL_Barrel", "ECAL_Endcap", "HCal_Endcap", "HCalBarrel", "HCalExtBarrel", "ECalFwd", "HCalFwd"}, "names of the calorimeters"};
+
+  /// Cell noise RMS (=1Sigma threshold) in GeV
   /// default values estimates for the LAr and TileCal from ATLAS Barrel: ATL-LARG-PUB-2008_002
   /// effective seed thresholds in topo-clustering of 7.5 and 11.5MeV in LAr and TileCal
-  /// EM Barrel Calorimeter
-  Gaudi::Property<double> m_ECalBThreshold{this, "ecalBarrelThreshold", 0.0075/4., "Cell threshold of EM Barrel Calorimeter in GeV"};
-  /// EM Endcap Calorimeter
-  Gaudi::Property<double> m_ECalECThreshold{this, "ecalEndcapThreshold", 0.0075/4., "Cell threshold of EM Endcap Calorimeter in GeV"};
-  /// EM Forward Calorimeter
-  Gaudi::Property<double> m_ECalFwdThreshold{this, "ecalFwdThreshold", 0.0075/4., "Cell threshold of EM Forward Calorimeter in GeV"};
-  /// Hadron Barrel Calorimeter
-  Gaudi::Property<double> m_HCalBThreshold{this, "hcalBarrelThreshold", 0.0115/4., "Cell threshold of Hadron Barrel Calorimeter in GeV"};
-  /// Hadron Endcap Calorimeter
-  Gaudi::Property<double> m_HCalECThreshold{this, "hcalEndcapThreshold", 0.0075/4., "Cell threshold of Hadron Endcap Calorimeter in GeV"};
-  /// Hadron Forward Calorimeter
-  Gaudi::Property<double> m_HCalFwdThreshold{this, "hcalFwdThreshold", 0.0075/4., "Cell threshold of Hadron Forward Calorimeter in GeV"};
+  Gaudi::Property<std::vector<double>> m_detectorsNoiseRMS{this, "detectorsNoiseRMS", {0.0075/4, 0.0075/4, 0.0115/4, 0.0115/4, 0.0115/4, 0.0075/4, 0.0075/4}, "Cell noise RMS in GeV"}; 
+
+  /// Set noise offset or not (false by default)
+  Gaudi::Property<bool> m_setNoiseOffset{this, "setNoiseOffset", false, "Set a noise offset per cell"};
+
+  /// Cell noise offset in GeV. 0 by default
+  Gaudi::Property<std::vector<double>> m_detectorsNoiseOffset{this, "detectorsNoiseOffset", {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}, "Cell noise offset in GeV"};
 
   /// Pointer to the geometry service
   SmartIF<IGeoSvc> m_geoSvc;
-  // Decoder
-  dd4hep::DDSegmentation::BitFieldCoder* m_decoder = new dd4hep::DDSegmentation::BitFieldCoder("system:4");
 
+  /// Decoder
+  dd4hep::DDSegmentation::BitFieldCoder* m_decoder = new dd4hep::DDSegmentation::BitFieldCoder("system:4");
 };
 
 #endif /* RECCALORIMETER_CONSTNOISETOOL_H */

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<int>(m_mu));
-      verbose() << " NUM CELLS = " << numCells << "   cluster noise const = " << getNoiseConstantPerCluster(newEta, numCells)
-                << " scaled to PU " << m_mu<< "  = " << getNoiseConstantPerCluster(newEta, numCells)* std::sqrt(static_cast<int>(m_mu)) << endmsg;
+      noise = getNoiseRMSPerCluster(newEta, numCells) * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
+      verbose() << " NUM CELLS = " << numCells << "   cluster noise const = " << getNoiseRMSPerCluster(newEta, numCells)
+                << " scaled to PU " << m_mu<< "  = " << getNoiseRMSPerCluster(newEta, numCells)* std::sqrt(static_cast<int>(m_mu)) << endmsg;
     } else {
       noise = m_constPileupNoise * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
     }
@@ -456,36 +456,36 @@ StatusCode CorrectECalBarrelSliWinCluster::initNoiseFromFile() {
   for (unsigned i = 0; i < 2; i++) {
     pileupParamHistoName = m_pileupHistoName + std::to_string(i);
     debug() << "Getting histogram with a name " << pileupParamHistoName << endmsg;
-    m_histoPileupConst.push_back(*dynamic_cast<TH1F*>(noiseFile->Get(pileupParamHistoName.c_str())));
-    if (m_histoPileupConst.at(i).GetNbinsX() < 1) {
+    m_histoPileupNoiseRMS.push_back(*dynamic_cast<TH1F*>(noiseFile->Get(pileupParamHistoName.c_str())));
+    if (m_histoPileupNoiseRMS.at(i).GetNbinsX() < 1) {
       error() << "Histogram  " << pileupParamHistoName
               << " has 0 bins! check the file with noise and the name of the histogram!" << endmsg;
       return StatusCode::FAILURE;
     }
   }
 
   // Check if we have same number of histograms (all layers) for pileup and electronics noise
-  if (m_histoPileupConst.size() == 0) {
+  if (m_histoPileupNoiseRMS.size() == 0) {
     error() << "No histograms with noise found!!!!" << endmsg;
     return StatusCode::FAILURE;
   }
   return StatusCode::SUCCESS;
 }
 
-double CorrectECalBarrelSliWinCluster::getNoiseConstantPerCluster(double aEta, uint aNumCells) {
+double CorrectECalBarrelSliWinCluster::getNoiseRMSPerCluster(double aEta, uint aNumCells) {
   double param0 = 0.;
   double param1 = 0.;
 
   // All histograms have same binning, all bins with same size
   // Using the histogram of the first parameter to get the bin size
   unsigned index = 0;
-  if (m_histoPileupConst.size() != 0) {
-    int Nbins = m_histoPileupConst.at(index).GetNbinsX();
+  if (m_histoPileupNoiseRMS.size() != 0) {
+    int Nbins = m_histoPileupNoiseRMS.at(index).GetNbinsX();
     double deltaEtaBin =
-        (m_histoPileupConst.at(index).GetBinLowEdge(Nbins) + m_histoPileupConst.at(index).GetBinWidth(Nbins) -
-         m_histoPileupConst.at(index).GetBinLowEdge(1)) /
+        (m_histoPileupNoiseRMS.at(index).GetBinLowEdge(Nbins) + m_histoPileupNoiseRMS.at(index).GetBinWidth(Nbins) -
+         m_histoPileupNoiseRMS.at(index).GetBinLowEdge(1)) /
         Nbins;
-    double etaFirtsBin = m_histoPileupConst.at(index).GetBinLowEdge(1);
+    double etaFirtsBin = m_histoPileupNoiseRMS.at(index).GetBinLowEdge(1);
     // find the eta bin for the cell
     int ibin = floor((fabs(aEta) - etaFirtsBin) / deltaEtaBin) + 1;
     verbose() << "Current eta = " << aEta << " bin = " << ibin << endmsg;
@@ -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;

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<edm4hep::ClusterCollection> 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<int> m_mu{this, "mu", 1000, "Pileup scenario"};
-  /// Histograms with pileup constants (index in array - parameter number: 0 or 1)
-  std::vector<TH1F> m_histoPileupConst;
+  /// Histograms with pileup RMS (index in array - parameter number: 0 or 1)
+  std::vector<TH1F> m_histoPileupNoiseRMS;
   /// Values of eta corresponding to the upstream correction parameters
   Gaudi::Property<std::vector<double>> m_etaValues{this, "etaValues", {0}, "Eta values"};
   /// Borders of the eta bins for the upstream correction (middle between eta values)

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<int>(m_mu));
-      verbose() << " NUM CELLS = " << numCells << "   cluster noise const = " << getNoiseConstantPerCluster(newEta, numCells)
-                << " scaled to PU " << m_mu<< "  = " << getNoiseConstantPerCluster(newEta, numCells)* std::sqrt(static_cast<int>(m_mu)) << endmsg;
+      noise = getNoiseRMSPerCluster(newEta, numCells) * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
+      verbose() << " NUM CELLS = " << numCells << "   cluster noise RMS = " << getNoiseRMSPerCluster(newEta, numCells)
+                << " scaled to PU " << m_mu<< "  = " << noise << endmsg;
     } else {
       noise = m_constPileupNoise * m_gauss.shoot() * std::sqrt(static_cast<int>(m_mu));
     }
@@ -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.;
 

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<edm4hep::ClusterCollection> m_inClusters{"clusters", Gaudi::DataHandle::Reader, this};
   /// Handle for corrected clusters (output collection)

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<uint64_t, double>& aCells) {
   std::for_each(aCells.begin(), aCells.end(),
-                [this](std::pair<const uint64_t, double>& p) { p.second += (m_gauss.shoot() * m_cellNoise); });
+                [this](std::pair<const uint64_t, double>& p) { p.second += (m_cellNoiseOffset + (m_gauss.shoot() * m_cellNoiseRMS)); });
 }
 
 void NoiseCaloCellsFlatTool::filterCellNoise(std::unordered_map<uint64_t, double>& aCells) {
   // Erase a cell if it has energy below a threshold
-  double threshold = m_filterThreshold * m_cellNoise;
+  double threshold = m_cellNoiseOffset + m_filterThreshold * m_cellNoiseRMS;
   auto it = aCells.begin();
   while ((it = std::find_if(it, aCells.end(), [&threshold](std::pair<const uint64_t, double>& p) {
             return bool(p.second < threshold);

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<uint64_t, double>& aCells) final;
-  /** @brief Remove cells with energy bellow threshold*sigma from the vector of cells
+  /** @brief Remove cells with energy below threshold*sigma from the vector of cells
    */
   virtual void filterCellNoise(std::unordered_map<uint64_t, double>& aCells) final;
 
 private:
-  /// Sigma of noise -- uniform noise per cell in GeV
-  Gaudi::Property<double> m_cellNoise{this, "cellNoise", 0.003, "uniform noise per cell in GeV"};
-  /// Energy threshold (Ecell < filterThreshold*m_cellNoise removed)
+  /// RMS of noise -- uniform RMS per cell in GeV
+  Gaudi::Property<double> m_cellNoiseRMS{this, "cellNoiseRMS", 0.003, "uniform noise RMS per cell in GeV"};
+  /// Offset of noise -- uniform offset per cell in GeV
+  Gaudi::Property<double> m_cellNoiseOffset{this, "cellNoiseOffset", 0.0, "uniform noise offset per cell in GeV"};
+  /// Energy threshold (Ecell < m_cellNoiseOffset + filterThreshold*m_cellNoiseRMS removed)
   Gaudi::Property<double> m_filterThreshold{
       this, "filterNoiseThreshold", 3,
-      "remove cells with energy bellow filterThreshold (threshold is multiplied by a cell noise sigma)"};
+      "remove cells with energy below offset + threshold * noise RMS"};
   /// Random Number Service
   IRndmGenSvc* m_randSvc;
   /// Gaussian random number generator used for smearing with a constant resolution (m_sigma)