CalorimeterHitDigi: add photon counter based readout

This implements modes to add additional signal smearing associated with
the scintillation process in calorimeters. The SiPM mode allows to study
effect of saturation for a given material, light collection efficiency
and type of sensor.

The part with binomial distribution sampling is a SiPM formalism from
Jin Huang's work done for sPhenix:
sPHENIX-Collaboration/coresoftware@f746f13

src/algorithms/calorimetry/CalorimeterHitDigi.cc

@@ -32,6 +32,7 @@
 #include <cstddef>
 #include <gsl/pointers>
 #include <limits>
+#include <map>
 #include <stdexcept>
 #include <string>
 #include <unordered_map>
@@ -122,6 +123,17 @@ void CalorimeterHitDigi::init() {
 
     auto& serviceSvc = algorithms::ServiceSvc::instance();
     corrMeanScale = serviceSvc.service<EvaluatorSvc>("EvaluatorSvc")->compile(m_cfg.corrMeanScale, hit_to_map);
+
+    std::map<std::string, readout_enum> readoutTypes {
+      {"simple", kSimpleReadout},
+      {"poisson_photon", kPoissonPhotonReadout},
+      {"sipm", kSipmReadout}
+    };
+    if (not readoutTypes.count(m_cfg.readoutType)) {
+        error("Invalid readoutType \"{}\"", m_cfg.readoutType);
+        throw std::runtime_error(fmt::format("Invalid readoutType \"{}\"", m_cfg.readoutType));
+    }
+    readoutType = readoutTypes.at(m_cfg.readoutType);
 }
 
 
@@ -205,14 +217,34 @@ void CalorimeterHitDigi::process(
                      std::pow(m_cfg.eRes[2] / (edep), 2)
                   )
                 : 0;
+
         double corrMeanScale_value = corrMeanScale(leading_hit);
 
         double ped = m_cfg.pedMeanADC + m_gaussian(m_generator) * m_cfg.pedSigmaADC;
 
         // Note: both adc and tdc values must be positive numbers to avoid integer wraparound
-        unsigned long long adc = std::max(std::llround(ped + edep * corrMeanScale_value * (1.0 + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC), 0LL);
+        unsigned long long adc;
         unsigned long long tdc = std::llround((time + m_gaussian(m_generator) * tRes) * stepTDC);
 
+        if (readoutType == kSimpleReadout) {
+            adc = std::max(std::llround(ped + edep * corrMeanScale_value * (1.0 + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC), 0LL);
+        } else if (readoutType == kPoissonPhotonReadout) {
+            const long long int n_photons_mean = edep * m_cfg.lightYield;
+            std::poisson_distribution<> n_photons_detected_dist(n_photons_mean);
+            const long long int n_photons_detected = n_photons_detected_dist(m_generator) * m_cfg.photonDetectionEfficiency;
+            const long long int n_max_photons = m_cfg.dyRangeADC * m_cfg.lightYield * m_cfg.photonDetectionEfficiency;
+            trace("n_photons_detected {}", n_photons_detected);
+            adc = std::max(std::llround(ped + n_photons_detected * corrMeanScale_value * (1.0 + eResRel) / n_max_photons * m_cfg.capADC), 0LL);
+        } else if (readoutType == kSipmReadout) {
+            const long long int n_photons = edep * m_cfg.lightYield;
+            std::binomial_distribution<> n_photons_detected_dist(n_photons, m_cfg.photonDetectionEfficiency);
+            const long long int n_photons_detected = n_photons_detected_dist(m_generator);
+            const long long int n_pixels_fired = m_cfg.numEffectiveSipmPixels * (1 - exp(- n_photons_detected / (double)m_cfg.numEffectiveSipmPixels));
+            const long long int n_max_photons = m_cfg.dyRangeADC * m_cfg.lightYield * m_cfg.photonDetectionEfficiency;
+            trace("n_photons_detected {}, n_pixels_fired {}, n_max_photons {}", n_photons_detected, n_pixels_fired, n_max_photons);
+            adc = std::max(std::llround(ped + n_pixels_fired * corrMeanScale_value * (1.0 + eResRel) / n_max_photons * m_cfg.capADC), 0LL);
+        }
+
         if (edep> 1.e-3) trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {} \t corrMeanScale: {}", edep, adc, time, m_cfg.capTime, tdc, corrMeanScale_value);
 
         rawhit.setCellID(leading_hit.getCellID());

src/algorithms/calorimetry/CalorimeterHitDigi.h

@@ -77,6 +77,9 @@ namespace eicrecon {
 
     dd4hep::IDDescriptor id_spec;
 
+    enum readout_enum { kSimpleReadout, kPoissonPhotonReadout, kSipmReadout };
+    enum readout_enum readoutType{kSimpleReadout};
+
   private:
     const algorithms::GeoSvc& m_geo = algorithms::GeoSvc::instance();
 

src/algorithms/calorimetry/CalorimeterHitDigiConfig.h

@@ -6,6 +6,9 @@
 #include <string>
 #include <vector>
 
+#include <Evaluator/DD4hepUnits.h>
+#include <edm4eic/unit_system.h>
+
 namespace eicrecon {
 
   struct CalorimeterHitDigiConfig {
@@ -16,6 +19,12 @@ namespace eicrecon {
     // single hit energy deposition threshold
     double                   threshold{1.0*dd4hep::keV};
 
+    // readout settings
+    std::string              readoutType{"simple"};
+    double                   lightYield{0. / edm4eic::unit::GeV};
+    double                   photonDetectionEfficiency; // (light collection efficiency) x (quantum efficiency)
+    unsigned long long       numEffectiveSipmPixels;
+
     // digitization settings
     unsigned int             capADC{1};
     double                   capTime{1000}; // dynamic range in ns

src/detectors/EEMC/EEMC.cc

@@ -43,6 +43,13 @@ extern "C" {
             .eRes = {0.0 * sqrt(dd4hep::GeV), 0.0, 0.0 * dd4hep::GeV},
             .tRes = 0.0 * dd4hep::ns,
             .threshold =  0.0 * dd4hep::MeV,  // Use ADC cut instead
+            .readoutType = "sipm",
+            .lightYield = 300. / dd4hep::MeV,
+            // See simulation study by A. Hoghmrtsyan https://indico.bnl.gov/event/20415/
+            // This includes quantum efficiency of the SiPM
+            .photonDetectionEfficiency = 17. / 300.,
+            // S14160-6015PS, 4 sensors per cell
+            .numEffectiveSipmPixels = 159565 * 4,
             .capADC = EcalEndcapN_capADC,
             .dyRangeADC = EcalEndcapN_dyRangeADC,
             .pedMeanADC = EcalEndcapN_pedMeanADC,

src/factories/calorimetry/CalorimeterHitDigi_factory.h

@@ -34,6 +34,10 @@ class CalorimeterHitDigi_factory : public JOmniFactory<CalorimeterHitDigi_factor
     ParameterRef<std::string> m_corrMeanScale {this, "scaleResponse", config().corrMeanScale};
     ParameterRef<std::vector<std::string>> m_fields {this, "signalSumFields", config().fields};
     ParameterRef<std::string> m_readout {this, "readoutClass", config().readout};
+    ParameterRef<std::string> m_readoutType {this, "readoutType", config().readoutType};
+    ParameterRef<double> m_lightYield {this, "lightYield", config().lightYield};
+    ParameterRef<double> m_photonDetectionEfficiency {this, "photonDetectionEfficiency", config().photonDetectionEfficiency};
+    ParameterRef<unsigned long long> m_numEffectiveSipmPixels {this, "numEffectiveSipmPixels", config().numEffectiveSipmPixels};
 
     Service<AlgorithmsInit_service> m_algorithmsInit {this};