Layer-theta sampling fraction integration for e/gamma

Detector/DetStudies/CMakeLists.txt

@@ -18,6 +18,7 @@ gaudi_add_module(DetStudies
                       EDM4HEP::edm4hep
                       ROOT::Core
                       ROOT::Hist
+                      ROOT::Tree
                 )
 
 install(TARGETS DetStudies

Detector/DetStudies/src/components/SamplingFractionInLayerThetaTree.cpp

@@ -0,0 +1,250 @@
+#include "SamplingFractionInLayerThetaTree.h"
+
+// FCCSW
+#include "k4Interface/IGeoSvc.h"
+
+// datamodel
+#include "edm4hep/SimCalorimeterHitCollection.h"
+
+#include "GaudiKernel/ITHistSvc.h"
+#include "TH2D.h"
+#include "TTree.h"
+
+// DD4hep
+#include "DD4hep/Detector.h"
+#include "DD4hep/Readout.h"
+
+#include <vector>
+
+DECLARE_COMPONENT(SamplingFractionInLayerThetaTree)
+
+SamplingFractionInLayerThetaTree::SamplingFractionInLayerThetaTree(const std::string& aName, ISvcLocator* aSvcLoc)
+    : Gaudi::Algorithm(aName, aSvcLoc), m_histSvc("THistSvc", aName), m_geoSvc("GeoSvc", aName), m_eventTree(nullptr),
+      m_summaryTree(nullptr), m_sfDistributionTree(nullptr), m_totalEnergyHist(nullptr), m_activeEnergyHist(nullptr),
+      m_eventCountHist(nullptr), m_samplingFractionHist(nullptr), m_sfDistribution(nullptr), m_event(0), m_layer(0),
+      m_theta(0), m_nHits(0), m_nActiveHits(0), m_summaryNHits(0), m_summaryNActiveHits(0), m_nEventBins(0),
+      m_totalEnergy(0.), m_activeEnergy(0.), m_samplingFraction(0.) {
+  declareProperty("deposits", m_deposits, "Energy deposits in sampling calorimeter (input)");
+}
+
+SamplingFractionInLayerThetaTree::~SamplingFractionInLayerThetaTree() {}
+
+StatusCode SamplingFractionInLayerThetaTree::initialize() {
+  if (Gaudi::Algorithm::initialize().isFailure()) {
+    return StatusCode::FAILURE;
+  }
+
+  if (!m_geoSvc) {
+    error() << "Unable to locate Geometry Service!" << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  if (m_geoSvc->getDetector()->readouts().find(m_readoutName) == m_geoSvc->getDetector()->readouts().end()) {
+    error() << "Readout <<" << m_readoutName << ">> does not exist." << endmsg;
+    return StatusCode::FAILURE;
+  }
+
+  const std::size_t nBins = static_cast<std::size_t>(m_numLayers) * static_cast<std::size_t>(m_numThetaBins);
+  m_sumEnergy.assign(nBins, 0.);
+  m_sumActiveEnergy.assign(nBins, 0.);
+  m_sumHits.assign(nBins, 0);
+  m_sumActiveHits.assign(nBins, 0);
+  m_sumEventBins.assign(nBins, 0);
+
+  if (m_writeEventTree) {
+    m_eventTree = new TTree("ecal_sf_layer_theta_tree", "Event-level layer-theta sampling fraction");
+    m_eventTree->Branch("event", &m_event, "event/I");
+    m_eventTree->Branch("layer", &m_layer, "layer/I");
+    m_eventTree->Branch("theta", &m_theta, "theta/I");
+    m_eventTree->Branch("total_energy", &m_totalEnergy, "total_energy/D");
+    m_eventTree->Branch("active_energy", &m_activeEnergy, "active_energy/D");
+    m_eventTree->Branch("sf", &m_samplingFraction, "sf/D");
+    m_eventTree->Branch("n_hits", &m_nHits, "n_hits/I");
+    m_eventTree->Branch("n_active_hits", &m_nActiveHits, "n_active_hits/I");
+
+    if (m_histSvc->regTree(m_treePath, m_eventTree).isFailure()) {
+      error() << "Couldn't register TTree at " << m_treePath << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+
+  if (m_writeSummaryTree) {
+    m_summaryTree = new TTree("ecal_sf_layer_theta_summary", "Final layer-theta sampling fraction summary");
+    m_summaryTree->Branch("layer", &m_layer, "layer/I");
+    m_summaryTree->Branch("theta", &m_theta, "theta/I");
+    m_summaryTree->Branch("total_energy", &m_totalEnergy, "total_energy/D");
+    m_summaryTree->Branch("active_energy", &m_activeEnergy, "active_energy/D");
+    m_summaryTree->Branch("sf", &m_samplingFraction, "sf/D");
+    m_summaryTree->Branch("n_hits", &m_summaryNHits, "n_hits/L");
+    m_summaryTree->Branch("n_active_hits", &m_summaryNActiveHits, "n_active_hits/L");
+    m_summaryTree->Branch("n_event_bins", &m_nEventBins, "n_event_bins/L");
+
+    if (m_histSvc->regTree(m_summaryTreePath, m_summaryTree).isFailure()) {
+      error() << "Couldn't register TTree at " << m_summaryTreePath << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+
+  if (m_writeMergeableHistograms) {
+    const double layerMin = static_cast<double>(m_firstLayerId) - 0.5;
+    const double layerMax = static_cast<double>(m_firstLayerId + static_cast<int>(m_numLayers)) - 0.5;
+    const double thetaMin = static_cast<double>(m_firstThetaId) - 0.5;
+    const double thetaMax = static_cast<double>(m_firstThetaId + static_cast<int>(m_numThetaBins)) - 0.5;
+
+    m_totalEnergyHist = new TH2D("ecal_sf_layer_theta_total_energy", "Total energy by layer and theta;layer;theta",
+                                 m_numLayers, layerMin, layerMax, m_numThetaBins, thetaMin, thetaMax);
+    m_activeEnergyHist = new TH2D("ecal_sf_layer_theta_active_energy", "Active energy by layer and theta;layer;theta",
+                                  m_numLayers, layerMin, layerMax, m_numThetaBins, thetaMin, thetaMax);
+    m_eventCountHist = new TH2D("ecal_sf_layer_theta_event_count", "Non-empty event bins by layer and theta;layer;theta",
+                                m_numLayers, layerMin, layerMax, m_numThetaBins, thetaMin, thetaMax);
+    m_samplingFractionHist = new TH2D("ecal_sf_layer_theta", "Final sampling fraction by layer and theta;layer;theta",
+                                      m_numLayers, layerMin, layerMax, m_numThetaBins, thetaMin, thetaMax);
+
+    if (m_histSvc->regHist(m_totalEnergyHistPath, m_totalEnergyHist).isFailure() ||
+        m_histSvc->regHist(m_activeEnergyHistPath, m_activeEnergyHist).isFailure() ||
+        m_histSvc->regHist(m_eventCountHistPath, m_eventCountHist).isFailure() ||
+        m_histSvc->regHist(m_samplingFractionHistPath, m_samplingFractionHist).isFailure()) {
+      error() << "Couldn't register layer-theta TH2D outputs" << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+
+  if (m_writeSfDistribution) {
+    if (m_sfDistributionMax <= m_sfDistributionMin) {
+      error() << "sfDistributionMax must be larger than sfDistributionMin." << endmsg;
+      return StatusCode::FAILURE;
+    }
+
+    const int bins[3] = {static_cast<int>(m_numLayers), static_cast<int>(m_numThetaBins),
+                         static_cast<int>(m_sfDistributionBins)};
+    const double min[3] = {static_cast<double>(m_firstLayerId) - 0.5, static_cast<double>(m_firstThetaId) - 0.5,
+                           m_sfDistributionMin};
+    const double max[3] = {static_cast<double>(m_firstLayerId + static_cast<int>(m_numLayers)) - 0.5,
+                           static_cast<double>(m_firstThetaId + static_cast<int>(m_numThetaBins)) - 0.5,
+                           m_sfDistributionMax};
+    m_sfDistribution = new THnSparseD("ecal_sf_layer_theta_sf_distribution",
+                                      "Event-level SF distribution by layer and theta;layer;theta;sf", 3, bins, min,
+                                      max);
+    m_sfDistributionTree = new TTree("ecal_sf_layer_theta_sf_distribution",
+                                     "THnSparseD event-level SF distribution by layer and theta");
+    m_sfDistributionTree->Branch("sf_distribution", &m_sfDistribution);
+    if (m_histSvc->regTree(m_sfDistributionPath, m_sfDistributionTree).isFailure()) {
+      error() << "Couldn't register THnSparseD holder TTree at " << m_sfDistributionPath << endmsg;
+      return StatusCode::FAILURE;
+    }
+  }
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode SamplingFractionInLayerThetaTree::execute(const EventContext& ctx) const {
+  auto decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder();
+
+  std::vector<std::vector<double>> sumE(m_numLayers, std::vector<double>(m_numThetaBins, 0.));
+  std::vector<std::vector<double>> sumEactive(m_numLayers, std::vector<double>(m_numThetaBins, 0.));
+  std::vector<std::vector<int>> nHits(m_numLayers, std::vector<int>(m_numThetaBins, 0));
+  std::vector<std::vector<int>> nActiveHits(m_numLayers, std::vector<int>(m_numThetaBins, 0));
+
+  const auto deposits = m_deposits.get();
+  for (const auto& hit : *deposits) {
+    dd4hep::DDSegmentation::CellID cID = hit.getCellID();
+    const int layerId = decoder->get(cID, m_layerFieldName);
+    const int thetaId = decoder->get(cID, m_thetaFieldName);
+    const int layerIndex = layerId - m_firstLayerId;
+    const int thetaIndex = thetaId - m_firstThetaId;
+
+    if (layerIndex < 0 || layerIndex >= static_cast<int>(m_numLayers) || thetaIndex < 0 ||
+        thetaIndex >= static_cast<int>(m_numThetaBins)) {
+      continue;
+    }
+
+    const double energy = hit.getEnergy();
+    sumE[layerIndex][thetaIndex] += energy;
+    ++nHits[layerIndex][thetaIndex];
+
+    const int activeField = decoder->get(cID, m_activeFieldName);
+    if (activeField == m_activeFieldValue) {
+      sumEactive[layerIndex][thetaIndex] += energy;
+      ++nActiveHits[layerIndex][thetaIndex];
+    }
+  }
+
+  m_event = static_cast<int>(ctx.evt());
+  for (unsigned int layerIndex = 0; layerIndex < m_numLayers; ++layerIndex) {
+    for (unsigned int thetaIndex = 0; thetaIndex < m_numThetaBins; ++thetaIndex) {
+      if (sumE[layerIndex][thetaIndex] <= 0.) {
+        continue;
+      }
+
+      m_layer = static_cast<int>(layerIndex) + m_firstLayerId;
+      m_theta = static_cast<int>(thetaIndex) + m_firstThetaId;
+      m_totalEnergy = sumE[layerIndex][thetaIndex];
+      m_activeEnergy = sumEactive[layerIndex][thetaIndex];
+      m_samplingFraction = m_activeEnergy / m_totalEnergy;
+      m_nHits = nHits[layerIndex][thetaIndex];
+      m_nActiveHits = nActiveHits[layerIndex][thetaIndex];
+
+      const std::size_t index = flatIndex(layerIndex, thetaIndex);
+      m_sumEnergy[index] += m_totalEnergy;
+      m_sumActiveEnergy[index] += m_activeEnergy;
+      m_sumHits[index] += m_nHits;
+      m_sumActiveHits[index] += m_nActiveHits;
+      ++m_sumEventBins[index];
+
+      if (m_writeMergeableHistograms) {
+        m_totalEnergyHist->Fill(m_layer, m_theta, m_totalEnergy);
+        m_activeEnergyHist->Fill(m_layer, m_theta, m_activeEnergy);
+        m_eventCountHist->Fill(m_layer, m_theta);
+      }
+
+      if (m_writeSfDistribution) {
+        const double coordinates[3] = {static_cast<double>(m_layer), static_cast<double>(m_theta), m_samplingFraction};
+        m_sfDistribution->Fill(coordinates);
+      }
+
+      if (m_writeEventTree) {
+        m_eventTree->Fill();
+      }
+    }
+  }
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode SamplingFractionInLayerThetaTree::finalize() {
+  for (unsigned int layerIndex = 0; layerIndex < m_numLayers; ++layerIndex) {
+    for (unsigned int thetaIndex = 0; thetaIndex < m_numThetaBins; ++thetaIndex) {
+      const std::size_t index = flatIndex(layerIndex, thetaIndex);
+      if (m_sumEnergy[index] <= 0.) {
+        continue;
+      }
+
+      m_layer = static_cast<int>(layerIndex) + m_firstLayerId;
+      m_theta = static_cast<int>(thetaIndex) + m_firstThetaId;
+      m_totalEnergy = m_sumEnergy[index];
+      m_activeEnergy = m_sumActiveEnergy[index];
+      m_samplingFraction = m_activeEnergy / m_totalEnergy;
+      m_summaryNHits = m_sumHits[index];
+      m_summaryNActiveHits = m_sumActiveHits[index];
+      m_nEventBins = m_sumEventBins[index];
+
+      if (m_writeSummaryTree) {
+        m_summaryTree->Fill();
+      }
+      if (m_writeMergeableHistograms) {
+        m_samplingFractionHist->SetBinContent(layerIndex + 1, thetaIndex + 1, m_samplingFraction);
+      }
+    }
+  }
+
+  if (m_writeSfDistribution) {
+    m_sfDistributionTree->Fill();
+  }
+
+  return Gaudi::Algorithm::finalize();
+}
+
+std::size_t SamplingFractionInLayerThetaTree::flatIndex(unsigned int layerIndex, unsigned int thetaIndex) const {
+  return static_cast<std::size_t>(layerIndex) * static_cast<std::size_t>(m_numThetaBins) +
+         static_cast<std::size_t>(thetaIndex);
+}

Detector/DetStudies/src/components/SamplingFractionInLayerThetaTree.h

@@ -0,0 +1,120 @@
+#ifndef DETSTUDIES_SAMPLINGFRACTIONINLAYERTHETATREE_H
+#define DETSTUDIES_SAMPLINGFRACTIONINLAYERTHETATREE_H
+
+// GAUDI
+#include "Gaudi/Algorithm.h"
+#include "GaudiKernel/ServiceHandle.h"
+
+// FCCSW
+#include "k4FWCore/DataHandle.h"
+
+#include "THnSparse.h"
+
+#include <cstddef>
+#include <vector>
+
+class IGeoSvc;
+class ITHistSvc;
+class TH2D;
+class TTree;
+
+// datamodel
+namespace edm4hep {
+class SimCalorimeterHitCollection;
+}
+
+/** @class SamplingFractionInLayerThetaTree SamplingFractionInLayerThetaTree.h
+ *
+ *  Accumulates layer-theta sampling fraction constants and writes final summary outputs.
+ *  It can also keep the legacy event-level TTree for workflows that still need per-event entries.
+ */
+class SamplingFractionInLayerThetaTree : public Gaudi::Algorithm {
+public:
+  explicit SamplingFractionInLayerThetaTree(const std::string&, ISvcLocator*);
+  virtual ~SamplingFractionInLayerThetaTree();
+
+  virtual StatusCode initialize() final;
+  virtual StatusCode execute(const EventContext&) const final;
+  virtual StatusCode finalize() final;
+
+private:
+  ServiceHandle<ITHistSvc> m_histSvc;
+  ServiceHandle<IGeoSvc> m_geoSvc;
+
+  mutable k4FWCore::DataHandle<edm4hep::SimCalorimeterHitCollection> m_deposits{"rec/caloHits",
+                                                                                Gaudi::DataHandle::Reader, this};
+
+  Gaudi::Property<std::string> m_activeFieldName{this, "activeFieldName", "", "Identifier of active material"};
+  Gaudi::Property<int> m_activeFieldValue{this, "activeFieldValue", 0, "Value of identifier for active material"};
+  Gaudi::Property<std::string> m_layerFieldName{this, "layerFieldName", "", "Identifier of layers"};
+  Gaudi::Property<std::string> m_thetaFieldName{this, "thetaFieldName", "theta", "Identifier of theta bins"};
+  Gaudi::Property<unsigned int> m_numLayers{this, "numLayers", 8, "Number of layers"};
+  Gaudi::Property<int> m_firstLayerId{this, "firstLayerId", 0, "First layer id included in output"};
+  Gaudi::Property<unsigned int> m_numThetaBins{this, "numThetaBins", 1024, "Number of theta bins"};
+  Gaudi::Property<int> m_firstThetaId{this, "firstThetaId", 0, "First theta id included in output"};
+  Gaudi::Property<std::string> m_readoutName{this, "readoutName", "", "Name of the detector readout"};
+  Gaudi::Property<bool> m_writeEventTree{this, "writeEventTree", true,
+                                         "Write the legacy event-level layer-theta TTree"};
+  Gaudi::Property<bool> m_writeSummaryTree{this, "writeSummaryTree", true,
+                                           "Write one final row per non-empty layer-theta bin"};
+  Gaudi::Property<bool> m_writeMergeableHistograms{this, "writeMergeableHistograms", true,
+                                                   "Write additive TH2D energy/count histograms for hadd merging"};
+  Gaudi::Property<bool> m_writeSfDistribution{this, "writeSfDistribution", false,
+                                              "Write THnSparseD(layer, theta, event-level SF) for per-bin fits"};
+  Gaudi::Property<unsigned int> m_sfDistributionBins{this, "sfDistributionBins", 200,
+                                                     "Number of bins in the THnSparseD SF axis"};
+  Gaudi::Property<double> m_sfDistributionMin{this, "sfDistributionMin", 0.,
+                                              "Minimum value of the THnSparseD SF axis"};
+  Gaudi::Property<double> m_sfDistributionMax{this, "sfDistributionMax", 1.0,
+                                              "Maximum value of the THnSparseD SF axis"};
+  Gaudi::Property<std::string> m_treePath{this, "treePath", "/rec/ecal_sf_layer_theta_tree",
+                                          "THistSvc path for the legacy event-level TTree"};
+  Gaudi::Property<std::string> m_summaryTreePath{this, "summaryTreePath", "/rec/ecal_sf_layer_theta_summary",
+                                                 "THistSvc path for the final summary TTree"};
+  Gaudi::Property<std::string> m_totalEnergyHistPath{this, "totalEnergyHistPath",
+                                                     "/rec/ecal_sf_layer_theta_total_energy",
+                                                     "THistSvc path for the additive total-energy TH2D"};
+  Gaudi::Property<std::string> m_activeEnergyHistPath{this, "activeEnergyHistPath",
+                                                      "/rec/ecal_sf_layer_theta_active_energy",
+                                                      "THistSvc path for the additive active-energy TH2D"};
+  Gaudi::Property<std::string> m_eventCountHistPath{this, "eventCountHistPath",
+                                                    "/rec/ecal_sf_layer_theta_event_count",
+                                                    "THistSvc path for the additive non-empty event-bin count TH2D"};
+  Gaudi::Property<std::string> m_samplingFractionHistPath{this, "samplingFractionHistPath",
+                                                          "/rec/ecal_sf_layer_theta",
+                                                          "THistSvc path for the final derived sampling-fraction TH2D"};
+  Gaudi::Property<std::string> m_sfDistributionPath{this, "sfDistributionPath",
+                                                   "/rec/ecal_sf_layer_theta_sf_distribution",
+                                                   "THistSvc path for the TTree holding the THnSparseD SF distribution"};
+
+  TTree* m_eventTree;
+  TTree* m_summaryTree;
+  TTree* m_sfDistributionTree;
+  TH2D* m_totalEnergyHist;
+  TH2D* m_activeEnergyHist;
+  TH2D* m_eventCountHist;
+  TH2D* m_samplingFractionHist;
+  THnSparseD* m_sfDistribution;
+
+  mutable int m_event;
+  mutable int m_layer;
+  mutable int m_theta;
+  mutable int m_nHits;
+  mutable int m_nActiveHits;
+  mutable long long m_summaryNHits;
+  mutable long long m_summaryNActiveHits;
+  mutable long long m_nEventBins;
+  mutable double m_totalEnergy;
+  mutable double m_activeEnergy;
+  mutable double m_samplingFraction;
+
+  mutable std::vector<double> m_sumEnergy;
+  mutable std::vector<double> m_sumActiveEnergy;
+  mutable std::vector<long long> m_sumHits;
+  mutable std::vector<long long> m_sumActiveHits;
+  mutable std::vector<long long> m_sumEventBins;
+
+  std::size_t flatIndex(unsigned int layerIndex, unsigned int thetaIndex) const;
+};
+
+#endif /* DETSTUDIES_SAMPLINGFRACTIONINLAYERTHETATREE_H */