From 4fcc8228a18e386d79f667a5e811bd3747d17ea9 Mon Sep 17 00:00:00 2001 From: Dmitry Kalinkin Date: Tue, 3 Oct 2023 17:41:01 -0400 Subject: [PATCH] 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: https://github.com/sPHENIX-Collaboration/coresoftware/commit/f746f130c11adfdb70b1978a5d3ebf1b258fe1e8 --- .../calorimetry/CalorimeterHitDigi.cc | 34 ++++++++++++++++++- .../calorimetry/CalorimeterHitDigi.h | 3 ++ .../calorimetry/CalorimeterHitDigiConfig.h | 9 +++++ src/detectors/EEMC/EEMC.cc | 7 ++++ .../calorimetry/CalorimeterHitDigi_factory.h | 4 +++ 5 files changed, 56 insertions(+), 1 deletion(-) diff --git a/src/algorithms/calorimetry/CalorimeterHitDigi.cc b/src/algorithms/calorimetry/CalorimeterHitDigi.cc index 95457082e9..d595b7aab1 100644 --- a/src/algorithms/calorimetry/CalorimeterHitDigi.cc +++ b/src/algorithms/calorimetry/CalorimeterHitDigi.cc @@ -32,6 +32,7 @@ #include #include #include +#include #include #include #include @@ -122,6 +123,17 @@ void CalorimeterHitDigi::init() { auto& serviceSvc = algorithms::ServiceSvc::instance(); corrMeanScale = serviceSvc.service("EvaluatorSvc")->compile(m_cfg.corrMeanScale, hit_to_map); + + std::map 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 * m_cfg.photonDetectionEfficiency; + std::poisson_distribution<> n_photons_detected_dist(n_photons_mean); + const long long int n_photons_detected = n_photons_detected_dist(m_generator); + 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()); diff --git a/src/algorithms/calorimetry/CalorimeterHitDigi.h b/src/algorithms/calorimetry/CalorimeterHitDigi.h index b2d919a55f..9f26021fa2 100644 --- a/src/algorithms/calorimetry/CalorimeterHitDigi.h +++ b/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(); diff --git a/src/algorithms/calorimetry/CalorimeterHitDigiConfig.h b/src/algorithms/calorimetry/CalorimeterHitDigiConfig.h index 2b25d55c46..f01decb7be 100644 --- a/src/algorithms/calorimetry/CalorimeterHitDigiConfig.h +++ b/src/algorithms/calorimetry/CalorimeterHitDigiConfig.h @@ -6,6 +6,9 @@ #include #include +#include +#include + 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 diff --git a/src/detectors/EEMC/EEMC.cc b/src/detectors/EEMC/EEMC.cc index 8908c7ac72..06550f9bc4 100644 --- a/src/detectors/EEMC/EEMC.cc +++ b/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, diff --git a/src/factories/calorimetry/CalorimeterHitDigi_factory.h b/src/factories/calorimetry/CalorimeterHitDigi_factory.h index 1fc28bc2da..4f69d93ac3 100644 --- a/src/factories/calorimetry/CalorimeterHitDigi_factory.h +++ b/src/factories/calorimetry/CalorimeterHitDigi_factory.h @@ -34,6 +34,10 @@ class CalorimeterHitDigi_factory : public JOmniFactory m_corrMeanScale {this, "scaleResponse", config().corrMeanScale}; ParameterRef> m_fields {this, "signalSumFields", config().fields}; ParameterRef m_readout {this, "readoutClass", config().readout}; + ParameterRef m_readoutType {this, "readoutType", config().readoutType}; + ParameterRef m_lightYield {this, "lightYield", config().lightYield}; + ParameterRef m_photonDetectionEfficiency {this, "photonDetectionEfficiency", config().photonDetectionEfficiency}; + ParameterRef m_numEffectiveSipmPixels {this, "numEffectiveSipmPixels", config().numEffectiveSipmPixels}; Service m_algorithmsInit {this};