From a87fcbee98fafc9b007ad4be8b99246c1ecbcfff Mon Sep 17 00:00:00 2001
From: "Blair D. Seidlitz" <bseidlitz@sphnxuser05.sdcc.bnl.gov>
Date: Thu, 4 Jun 2026 23:29:25 -0400
Subject: [PATCH] CEMC SiPM pixel occupancy correction waveform sim

---
 .../g4waveformsim/CaloWaveformSim.cc          | 55 ++++++++++++++++++-
 .../g4waveformsim/CaloWaveformSim.h           |  2 +
 2 files changed, 55 insertions(+), 2 deletions(-)

diff --git a/simulation/g4simulation/g4waveformsim/CaloWaveformSim.cc b/simulation/g4simulation/g4waveformsim/CaloWaveformSim.cc
index cb2ec42d8b..ce2a7626f0 100644
--- a/simulation/g4simulation/g4waveformsim/CaloWaveformSim.cc
+++ b/simulation/g4simulation/g4waveformsim/CaloWaveformSim.cc
@@ -36,6 +36,7 @@
 #include <TTree.h>
 #include <algorithm>
 #include <cassert>
+#include <cmath>
 #include <sstream>
 #include <string>
 
@@ -55,6 +56,11 @@ CaloWaveformSim::~CaloWaveformSim()
   gsl_rng_free(m_RandomGenerator);
 }
 
+void CaloWaveformSim::set_use_sipm_occupancy(bool use_sipm_occupancy)
+{
+  m_use_sipm_occupancy = use_sipm_occupancy;
+}
+
 int CaloWaveformSim::InitRun(PHCompositeNode *topNode)
 {
   // Initialize random generator
@@ -295,6 +301,7 @@ int CaloWaveformSim::process_event(PHCompositeNode *topNode)
   }
 
   std::map<unsigned int,float> tbt_smear;
+  std::map<unsigned int, double> tower_photon_count_mean;
 
 
   // loop over hits
@@ -314,9 +321,11 @@ int CaloWaveformSim::process_event(PHCompositeNode *topNode)
     float correction = 1.;
     maphitetaphi(hit, etabin, phibin, correction);
     unsigned int key = encode_tower(etabin, phibin);
+    unsigned int tower_index = decode_tower(key);
     float calibconst = cdbttree->GetFloatValue(key, m_fieldname);
     float e_vis = hit->get_light_yield();
     e_vis *= correction;
+
     if (m_smear_const)
     {
       auto it = tbt_smear.find(key);
@@ -327,10 +336,28 @@ int CaloWaveformSim::process_event(PHCompositeNode *topNode)
       }
       else
       {
-        tbt_smear[key] =  1.0+ gsl_ran_gaussian(m_RandomGenerator,factor_const);
+        float val =  1.0+ gsl_ran_gaussian(m_RandomGenerator,factor_const);
+        if(val < 0.0f)
+        {
+          val = 0;
+        }
+        tbt_smear[key] = val; 
         e_vis *= tbt_smear[key];
       }
     }
+
+    if (m_use_sipm_occupancy && m_dettype == CaloTowerDefs::CEMC)
+    {
+      constexpr double kSamplingFraction = 2e-2;
+      constexpr double kPhotoelectronsPerGeV = 500.;
+      constexpr double kPhotonElecYieldVisibleGeV = kPhotoelectronsPerGeV / kSamplingFraction;
+      const double photon_count_mean = static_cast<double>(e_vis) * kPhotonElecYieldVisibleGeV;
+      if (photon_count_mean > 0.)
+      {
+        tower_photon_count_mean[tower_index] += photon_count_mean;
+      }
+    }
+
     float e_dep = e_vis / m_sampling_fraction;
     float ADC = (calibconst != 0) ? e_dep / calibconst : 0.;
     ADC *= m_gain;
@@ -351,7 +378,6 @@ int CaloWaveformSim::process_event(PHCompositeNode *topNode)
     }
 
     float t0 = hit->get_t(0) / m_sampletime;
-    unsigned int tower_index = decode_tower(key);
     // here I will add the truth matching part
     //  for the cell reco, the truth matching info relys on edep not light yield, I will be consistent here :)
     TowerInfo *tower = m_CaloWaveformContainer->get_tower_at_channel(tower_index);
@@ -371,6 +397,31 @@ int CaloWaveformSim::process_event(PHCompositeNode *topNode)
     }
   }
 
+  if (m_use_sipm_occupancy && m_dettype == CaloTowerDefs::CEMC)
+  {
+    constexpr double kSiPMEffectivePixel = 40000 * 4.;
+    for (const auto &entry : tower_photon_count_mean)
+    {
+      const unsigned int tower_index = entry.first;
+      const double photon_count_mean = entry.second;
+      if (photon_count_mean <= 0. || tower_index >= m_waveforms.size())
+      {
+        continue;
+      }
+
+      const double poisson_param_per_pixel = photon_count_mean / kSiPMEffectivePixel;
+      const double expected_active_pixels =
+          kSiPMEffectivePixel * (1. - std::exp(-poisson_param_per_pixel));
+      const double occupancy_ratio =
+          std::max(0., std::min(1., expected_active_pixels / photon_count_mean));
+      //std::cout << "occupancy_ratio: " << occupancy_ratio << std::endl;
+      for (int isample = 0; isample < m_nsamples; ++isample)
+      {
+        m_waveforms.at(tower_index).at(isample) *= occupancy_ratio;
+      }
+    }
+  }
+
   // do noise here and add to waveform
 
   if (m_noiseType == NoiseType::NOISE_TREE)
diff --git a/simulation/g4simulation/g4waveformsim/CaloWaveformSim.h b/simulation/g4simulation/g4waveformsim/CaloWaveformSim.h
index 80290a5321..573b8ce5d3 100644
--- a/simulation/g4simulation/g4waveformsim/CaloWaveformSim.h
+++ b/simulation/g4simulation/g4waveformsim/CaloWaveformSim.h
@@ -153,6 +153,7 @@ class CaloWaveformSim : public SubsysReco
 
   // Light collection model access
   LightCollectionModel &get_light_collection_model() { return light_collection_model; }
+  void set_use_sipm_occupancy(bool use_sipm_occupancy);
 
  private:
   CaloTowerDefs::DetectorSystem m_dettype{CaloTowerDefs::CEMC};
@@ -227,6 +228,7 @@ class CaloWaveformSim : public SubsysReco
   CDBTTree *cdbttree_time{nullptr}, *cdbttree_MC_time{nullptr};
   TProfile *h_template{nullptr};
   LightCollectionModel light_collection_model;
+  bool m_use_sipm_occupancy{true};
 
   NoiseType m_noiseType{NOISE_TREE};