Skip to content

CEMC SiPM pixel occupancy correction waveform sim #4286

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
bseidlit wants to merge 2 commits into
base: master
Choose a base branch
from
+55 −2
Draft

CEMC SiPM pixel occupancy correction waveform sim #4286

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a87fcbe
CEMC SiPM pixel occupancy correction waveform sim
Jun 5, 2026
01b0d45
Merge branch 'sPHENIX-Collaboration:master' into g4waveform
bseidlit Jun 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
55 changes: 53 additions & 2 deletions simulation/g4simulation/g4waveformsim/CaloWaveformSim.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
#include <TTree.h>
#include <algorithm>
#include <cassert>
#include <cmath>
#include <sstream>
#include <string>

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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);
Expand All @@ -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;
}
}
Comment on lines +349 to +359

@coderabbitai coderabbitai Bot Jun 5, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major | ⚡ Quick win

Hardcoded kSamplingFraction may diverge from configurable m_sampling_fraction.

The photon count calculation uses a hardcoded kSamplingFraction = 2e-2, but the energy deposit calculation at line 361 uses the configurable m_sampling_fraction member. If a user overrides the sampling fraction via set_sampling_fraction(), the photon statistics and energy calculations will be inconsistent.

Similarly, kPhotoelectronsPerGeV (500) has no corresponding configurable parameter, making it difficult to tune this model for different detector configurations.

Suggested fix: use member variable 
     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 kPhotonElecYieldVisibleGeV = kPhotoelectronsPerGeV / m_sampling_fraction;
       const double photon_count_mean = static_cast<double>(e_vis) * kPhotonElecYieldVisibleGeV;
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
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;
}
}
if (m_use_sipm_occupancy && m_dettype == CaloTowerDefs::CEMC)
{
constexpr double kPhotoelectronsPerGeV = 500.;
const double kPhotonElecYieldVisibleGeV = kPhotoelectronsPerGeV / m_sampling_fraction;
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;
Expand All @@ -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);
Expand All @@ -371,6 +397,31 @@ int CaloWaveformSim::process_event(PHCompositeNode *topNode)
}
}

if (m_use_sipm_occupancy && m_dettype == CaloTowerDefs::CEMC)
{
constexpr double kSiPMEffectivePixel = 40000 * 4.;

@pinkenburg pinkenburg Jun 6, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

there is already a definition in RawTowerDigitizer.h

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)
Expand Down
2 changes: 2 additions & 0 deletions simulation/g4simulation/g4waveformsim/CaloWaveformSim.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -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};

@coderabbitai coderabbitai Bot Jun 5, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⚠️ Potential issue | 🟠 Major | ⚡ Quick win

Default true silently changes CEMC simulation behavior.

The new SiPM occupancy scaling feature defaults to enabled, which will change waveform amplitudes for all existing CEMC simulations without user intervention. Per collaboration guidelines, simulation behavior changes should document expected analysis impact and whether reprocessing is required.

Consider defaulting to false for backwards compatibility, requiring explicit opt-in:

Suggested change 
-  bool m_use_sipm_occupancy{true};
+  bool m_use_sipm_occupancy{false};

Alternatively, document the expected impact on existing physics analyses in the PR description. Based on learnings: "If the PR changes simulation behavior, ensure the description states expected analysis impact and whether reprocessing is required."

📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
bool m_use_sipm_occupancy{true};
bool m_use_sipm_occupancy{false};


NoiseType m_noiseType{NOISE_TREE};

Expand Down