Skip to content

sEPD Q Vector Calibration (Generation): Charge Clamping and Ring 0 Exclusion #4245

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

Merged
pinkenburg merged 4 commits into from
Apr 12, 2026
Merged

sEPD Q Vector Calibration (Generation): Charge Clamping and Ring 0 Exclusion #4245

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
eac9a47
sEPD TreeGen: Remove hardcoded sEPD input node name
Steepspace Mar 31, 2026
316fbf2
QVecCalib - Add EP Resolution QA Hist
Steepspace Mar 31, 2026
19956d8
QVecCalib: sEPD Q-vector Calibration and Channel Selection
Steepspace Apr 11, 2026
169c06b
QVecCalib: Code Review Fixes
Steepspace Apr 11, 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
218 changes: 36 additions & 182 deletions calibrations/sepd/sepd_eventplanecalib/QVecCalib.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -111,13 +111,6 @@ int QVecCalib::process_QA_hist()
return ret;
}

// Get List of Bad Channels
ret = process_bad_channels(file);
if (ret)
{
return ret;
}

// cleanup
file->Close();
delete file;
Expand Down Expand Up @@ -193,138 +186,6 @@ int QVecCalib::process_sEPD_event_thresholds(TFile* file)
return Fun4AllReturnCodes::EVENT_OK;
}

int QVecCalib::process_bad_channels(TFile* file)
{
Fun4AllServer *se = Fun4AllServer::instance();

std::string sepd_charge_hist = "hSEPD_Charge";

TProfile *hSEPD_Charge{nullptr};
file->GetObject(sepd_charge_hist.c_str(), hSEPD_Charge);

// Check if the hist is stored in the file
if (hSEPD_Charge == nullptr)
{
std::cout << PHWHERE << "Error! Cannot find hist: " << sepd_charge_hist << ", in file: " << file->GetName() << std::endl;
return Fun4AllReturnCodes::ABORTRUN;
}

int rbins = 16;
int bins_charge = 40;

h2SEPD_South_Charge_rbin = new TH2F("h2SEPD_South_Charge_rbin",
"sEPD South; r_{bin}; Avg Charge",
rbins, -0.5, rbins - 0.5,
bins_charge, 0, bins_charge);

h2SEPD_North_Charge_rbin = new TH2F("h2SEPD_North_Charge_rbin",
"sEPD North; r_{bin}; Avg Charge",
rbins, -0.5, rbins - 0.5,
bins_charge, 0, bins_charge);

h2SEPD_South_Charge_rbinv2 = new TH2F("h2SEPD_South_Charge_rbinv2",
"sEPD South; r_{bin}; Avg Charge",
rbins, -0.5, rbins - 0.5,
bins_charge, 0, bins_charge);

h2SEPD_North_Charge_rbinv2 = new TH2F("h2SEPD_North_Charge_rbinv2",
"sEPD North; r_{bin}; Avg Charge",
rbins, -0.5, rbins - 0.5,
bins_charge, 0, bins_charge);

hSEPD_Bad_Channels = new TProfile("h_sEPD_Bad_Channels", "sEPD Bad Channels; Channel; Status", QVecShared::SEPD_CHANNELS, -0.5, QVecShared::SEPD_CHANNELS-0.5);

se->registerHisto(h2SEPD_South_Charge_rbin);
se->registerHisto(h2SEPD_North_Charge_rbin);
se->registerHisto(h2SEPD_South_Charge_rbinv2);
se->registerHisto(h2SEPD_North_Charge_rbinv2);
se->registerHisto(hSEPD_Bad_Channels);

for (int channel = 0; channel < QVecShared::SEPD_CHANNELS; ++channel)
{
unsigned int key = TowerInfoDefs::encode_epd(channel);
int rbin = TowerInfoDefs::get_epd_rbin(key);
unsigned int arm = TowerInfoDefs::get_epd_arm(key);

double avg_charge = hSEPD_Charge->GetBinContent(channel + 1);

auto* h2 = (arm == 0) ? h2SEPD_South_Charge_rbin : h2SEPD_North_Charge_rbin;

h2->Fill(rbin, avg_charge);
}

auto* hSpx = h2SEPD_South_Charge_rbin->ProfileX("hSpx", 2, -1, "s");
auto* hNpx = h2SEPD_North_Charge_rbin->ProfileX("hNpx", 2, -1, "s");

int ctr_dead = 0;
int ctr_hot = 0;
int ctr_cold = 0;

for (int channel = 0; channel < QVecShared::SEPD_CHANNELS; ++channel)
{
unsigned int key = TowerInfoDefs::encode_epd(channel);
int rbin = TowerInfoDefs::get_epd_rbin(key);
unsigned int arm = TowerInfoDefs::get_epd_arm(key);

auto* h2 = (arm == 0) ? h2SEPD_South_Charge_rbinv2 : h2SEPD_North_Charge_rbinv2;
auto* hprof = (arm == 0) ? hSpx : hNpx;

double charge = hSEPD_Charge->GetBinContent(channel + 1);
double mean_charge = hprof->GetBinContent(rbin + 1);
double sigma = hprof->GetBinError(rbin + 1);
double zscore = 0.0;

if (sigma > 0)
{
zscore = (charge - mean_charge) / sigma;
}

if (charge < m_sEPD_min_avg_charge_threshold || std::abs(zscore) > m_sEPD_sigma_threshold)
{
m_bad_channels.insert(channel);

std::string type;
QVecShared::ChannelStatus status_fill;

// dead channel
if (charge == 0)
{
type = "Dead";
status_fill = QVecShared::ChannelStatus::Dead;
++ctr_dead;
}
// hot channel
else if (zscore > m_sEPD_sigma_threshold)
{
type = "Hot";
status_fill = QVecShared::ChannelStatus::Hot;
++ctr_hot;
}
// cold channel
else
{
type = "Cold";
status_fill = QVecShared::ChannelStatus::Cold;
++ctr_cold;
}

hSEPD_Bad_Channels->Fill(channel, static_cast<int>(status_fill));
std::cout << std::format("{:4} Channel: {:3d}, arm: {}, rbin: {:2d}, Mean: {:5.2f}, Charge: {:5.2f}, Z-Score: {:5.2f}",
type, channel, arm, rbin, mean_charge, charge, zscore) << std::endl;
}
else
{
h2->Fill(rbin, charge);
}
}

std::cout << "Total Bad Channels: " << m_bad_channels.size() << ", Dead: "
<< ctr_dead << ", Hot: " << ctr_hot << ", Cold: " << ctr_cold << std::endl;

std::cout << "Finished processing Hot sEPD channels" << std::endl;
return Fun4AllReturnCodes::EVENT_OK;
}

void QVecCalib::init_hists()
{
unsigned int bins_psi = 126;
Expand Down Expand Up @@ -361,6 +222,14 @@ void QVecCalib::init_hists()
m_hists2D[name_N] = new TH2F(name_N.c_str(), title_N.c_str(), m_cent_bins, m_cent_low, m_cent_high, bins_psi, psi_low, psi_high);
m_hists2D[name_NS] = new TH2F(name_NS.c_str(), title_NS.c_str(), m_cent_bins, m_cent_low, m_cent_high, bins_psi, psi_low, psi_high);

if (m_pass == Pass::ApplyFlattening)
{
std::string name_EP_res = std::format("hEP_res_{}", n);
std::string title_EP_res = std::format("; Centrality [%]; #LTRe(Q^{{S}}_{{{0}}} Q^{{N*}}_{{{0}}}) / (|Q^{{S}}_{{{0}}}||Q^{{N}}_{{{0}}}|)#GT", n);

m_profiles[name_EP_res] = new TProfile(name_EP_res.c_str(), title_EP_res.c_str(), m_cent_bins, m_cent_low, m_cent_high);
}

// South, North
for (auto det : m_subdetectors)
{
Expand Down Expand Up @@ -723,6 +592,8 @@ void QVecCalib::prepare_flattening_hists()
std::string psi_N_name = std::format("h2_sEPD_Psi_N_{}_corr2", n);
std::string psi_NS_name = std::format("h2_sEPD_Psi_NS_{}_corr2", n);

std::string EP_res_name = std::format("hEP_res_{}", n);

FlatteningHists h;

h.S_x_corr2_avg = m_profiles.at(S_x_corr2_avg_name);
Expand All @@ -746,6 +617,8 @@ void QVecCalib::prepare_flattening_hists()
h.Psi_N_corr2 = m_hists2D.at(psi_N_name);
h.Psi_NS_corr2 = m_hists2D.at(psi_NS_name);

h.EP_res = m_profiles.at(EP_res_name);

se->registerHisto(h.S_x_corr2_avg);
se->registerHisto(h.S_y_corr2_avg);
se->registerHisto(h.N_x_corr2_avg);
Expand All @@ -767,6 +640,8 @@ void QVecCalib::prepare_flattening_hists()
se->registerHisto(h.Psi_N_corr2);
se->registerHisto(h.Psi_NS_corr2);

se->registerHisto(h.EP_res);

m_flattening_hists.push_back(h);
}
}
Expand Down Expand Up @@ -908,6 +783,11 @@ void QVecCalib::process_flattening(double cent, size_t h_idx, const QVecShared::
double psi_N = std::atan2(q_N_corr2.y, q_N_corr2.x);
double psi_NS = std::atan2(q_NS_corr2.y, q_NS_corr2.x);

double SP_QS_QN = q_S_corr2.x * q_N_corr2.x + q_S_corr2.y * q_N_corr2.y;
double norm_S = std::sqrt(q_S_corr2.x * q_S_corr2.x + q_S_corr2.y * q_S_corr2.y);
double norm_N = std::sqrt(q_N_corr2.x * q_N_corr2.x + q_N_corr2.y * q_N_corr2.y);
double EP_res = (norm_S && norm_N) ? SP_QS_QN / (norm_S * norm_N) : 0;

h.S_x_corr2_avg->Fill(cent, q_S_corr2.x);
h.S_y_corr2_avg->Fill(cent, q_S_corr2.y);
h.N_x_corr2_avg->Fill(cent, q_N_corr2.x);
Expand All @@ -927,6 +807,8 @@ void QVecCalib::process_flattening(double cent, size_t h_idx, const QVecShared::
h.Psi_S_corr2->Fill(cent, psi_S);
h.Psi_N_corr2->Fill(cent, psi_N);
h.Psi_NS_corr2->Fill(cent, psi_NS);

h.EP_res->Fill(cent, EP_res);
}

bool QVecCalib::process_sEPD()
Expand All @@ -939,14 +821,27 @@ bool QVecCalib::process_sEPD()
{
double charge = m_evtdata->get_sepd_charge(channel);

// Skip Bad Channels
if (m_bad_channels.contains(channel) || charge <= 0)
// Skip Noise
if (charge <= m_sEPD_noise_threshold)
{
continue;
}

// Clamp on high charge threshold
if (m_sEPD_charge_threshold > 0 && charge > m_sEPD_charge_threshold)
{
charge = m_sEPD_charge_threshold;
}

unsigned int key = TowerInfoDefs::encode_epd(channel);
unsigned int arm = TowerInfoDefs::get_epd_arm(key);
int rbin = TowerInfoDefs::get_epd_rbin(key);

// Skip Innermost Ring
if (rbin == 0)
{
continue;
}
Comment on lines 836 to +844
Copy link

Copilot AI Apr 11, 2026

Choose a reason for hiding this comment

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

Ring 0 is excluded from the Q-vector accumulation here, but Ring 0 still contributes to EventPlaneData::sepd_totalcharge (set upstream) which is used by process_event_check() and the QA-derived charge thresholds. If the intent is to mitigate Ring 0 beam-background effects, the event-selection charge (and/or the QA histogram used to derive its bounds) should be computed with the same Ring 0 exclusion (and clamping/noise definition) as the Q-vector loop.

Copilot uses AI. Check for mistakes.

// arm = 0: South
// arm = 1: North
Expand Down Expand Up @@ -1273,50 +1168,9 @@ void QVecCalib::write_cdb()
std::cout << "Info: Directory " << m_cdb_output_dir << " already exists." << std::endl;
}

write_cdb_BadTowers();
write_cdb_EventPlane();
}

void QVecCalib::write_cdb_BadTowers()
{
std::cout << "Writing Bad Towers CDB" << std::endl;

std::string payload = "SEPD_HotMap";
std::string fieldname_status = "status";
std::string fieldname_sigma = "SEPD_sigma";
std::string output_file = std::format("{}/{}-{}-{}.root", m_cdb_output_dir, payload, m_dst_tag, m_runnumber);

CDBTTree cdbttree(output_file);

for (int channel = 0; channel < QVecShared::SEPD_CHANNELS; ++channel)
{
unsigned int key = TowerInfoDefs::encode_epd(channel);
int status = hSEPD_Bad_Channels->GetBinContent(channel+1);

float sigma = 0;

// Hot
if (status == static_cast<int>(QVecShared::ChannelStatus::Hot))
{
sigma = SIGMA_HOT;
}

// Cold
else if (status == static_cast<int>(QVecShared::ChannelStatus::Cold))
{
sigma = SIGMA_COLD;
}

cdbttree.SetIntValue(key, fieldname_status, status);
cdbttree.SetFloatValue(key, fieldname_sigma, sigma);
}

std::cout << "Saving CDB: " << payload << " to " << output_file << std::endl;

cdbttree.Commit();
cdbttree.WriteCDBTTree();
}

void QVecCalib::write_cdb_EventPlane()
{
std::cout << "Writing Event Plane CDB" << std::endl;
Expand Down
44 changes: 16 additions & 28 deletions calibrations/sepd/sepd_eventplanecalib/QVecCalib.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@

#include <fun4all/SubsysReco.h>

#include <algorithm>
#include <array>
#include <map>
#include <stdexcept>
Expand Down Expand Up @@ -96,6 +97,16 @@ class QVecCalib : public SubsysReco
m_cdb_output_dir = cdb_dir;
}

void set_charge_threshold(double threshold)
{
m_sEPD_charge_threshold = std::max(0.0, threshold);
}
Comment thread
Steepspace marked this conversation as resolved.

void set_noise_threshold(double threshold)
{
m_sEPD_noise_threshold = threshold;
Comment thread
Steepspace marked this conversation as resolved.
Comment thread
Steepspace marked this conversation as resolved.
}

private:
static Pass validate_pass(int pass)
{
Expand Down Expand Up @@ -222,34 +233,27 @@ class QVecCalib : public SubsysReco
TProfile* NS_yy_corr_avg{nullptr};
TProfile* NS_xy_corr_avg{nullptr};

TProfile* EP_res{nullptr};

TH2* Psi_S_corr2{nullptr};
TH2* Psi_N_corr2{nullptr};
TH2* Psi_NS_corr2{nullptr};
};

// sEPD Bad Channels
std::unordered_set<int> m_bad_channels;

double m_sEPD_min_avg_charge_threshold{1};
double m_sEPD_sigma_threshold{3};

double m_sEPD_charge_threshold{50};
double m_sEPD_noise_threshold{0.5};

// Hists
TH1* hCentrality{nullptr};

TH2* h2SEPD_Charge{nullptr};
TH2* h2SEPD_Chargev2{nullptr};

TH2* h2SEPD_South_Charge_rbin{nullptr};
TH2* h2SEPD_North_Charge_rbin{nullptr};

TH2* h2SEPD_South_Charge_rbinv2{nullptr};
TH2* h2SEPD_North_Charge_rbinv2{nullptr};

TProfile* hSEPD_Charge_Min{nullptr};
TProfile* hSEPD_Charge_Max{nullptr};

TProfile* hSEPD_Bad_Channels{nullptr};

std::map<std::string, TH2*> m_hists2D;
std::map<std::string, TProfile*> m_profiles;

Expand Down Expand Up @@ -397,14 +401,6 @@ class QVecCalib : public SubsysReco
*/
int process_QA_hist();

/**
* @brief Identifies and catalogs "Bad" (Hot, Cold, or Dead) sEPD channels.
* * Uses a reference charge histogram to compute Z-scores based on mean charge
* per radial bin. Channels exceeding the sigma threshold are added to the internal exclusion set.
* * @param file Pointer to the open TFile containing QA histograms.
*/
int process_bad_channels(TFile* file);

/**
* @brief Establishes sEPD charge-cut thresholds for event selection.
* * Uses the 2D total charge vs. centrality distribution to derive mean and
Expand All @@ -422,14 +418,6 @@ class QVecCalib : public SubsysReco
* * @param output_dir The filesystem directory where the .root payload will be saved.
*/
void write_cdb_EventPlane();

/**
* @brief Writes the Hot/Cold tower status map to a CDB-formatted TTree.
* * Encodes sEPD channel indices into TowerInfo keys and maps status codes (1=Dead,
* 2=Hot, 3=Cold) to the final database payload.
* * @param output_dir The filesystem directory where the .root payload will be saved.
*/
void write_cdb_BadTowers();
};

#endif // SEPDEVENTPLANECALIB_QVECCALIB_H
Loading