Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
77 changes: 64 additions & 13 deletions src/programs/Simulation/mcsmear/FDCSmearer.cc
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,14 @@

#include <sstream>

bool fdc_config_t::FDC_EFFVSDOCA(1);
double fdc_config_t::FDC_CATHODE_SIGMA_SET(-1.);
double fdc_config_t::FDC_DRIFT_SMEAR_NS(0.);
double fdc_config_t::FDC_GAIN_SCALE(1.);
double fdc_config_t::FDC_GAIN_PKG[4] = {1., 1., 1., 1.};
double fdc_config_t::FDC_DOCA_ATTEN(0.);
double fdc_config_t::FDC_GAIN_SIGMA(0.);
double fdc_config_t::FDC_STRIP_THR_SET(-1.);
double fdc_config_t::FDC_WIRE_THR_KEV_SET(-1.);

//-----------
// fdc_config_t (constructor)
Expand Down Expand Up @@ -40,7 +47,12 @@ fdc_config_t::fdc_config_t(const std::shared_ptr<const JEvent>& event)
// Calculate ped noise level based on position resolution
// FDC_PED_NOISE = -0.004594 + 0.008711*FDC_CATHODE_SIGMA +
// 0.000010*FDC_CATHODE_SIGMA*FDC_CATHODE_SIGMA; //pC
if (FDC_CATHODE_SIGMA_SET >= 0) FDC_CATHODE_SIGMA = FDC_CATHODE_SIGMA_SET;
FDC_PED_NOISE = -0.0938 + 0.0485*FDC_CATHODE_SIGMA;
jout << "FDC_DRIFT_SMEAR_NS = " << FDC_DRIFT_SMEAR_NS << endl;
jout << "FDC_GAIN_SCALE = " << FDC_GAIN_SCALE << " x package " << FDC_GAIN_PKG[0] << "/" << FDC_GAIN_PKG[1] << "/" << FDC_GAIN_PKG[2] << "/" << FDC_GAIN_PKG[3] << ", FDC_DOCA_ATTEN = " << FDC_DOCA_ATTEN << ", FDC_GAIN_SIGMA = " << FDC_GAIN_SIGMA << ", FDC_STRIP_THRESHOLD = " << FDC_STRIP_THR_SET << " pC, FDC_WIRE_THRESHOLD_KEV = " << FDC_WIRE_THR_KEV_SET << " (< 0: CCDB)" << endl;
jout << "FDC_THRESH_KEV = " << FDC_THRESH_KEV << " (anode wire threshold)" << endl;
jout << "FDC_CATHODE_SIGMA = " << FDC_CATHODE_SIGMA << " um, FDC_PED_NOISE = " << FDC_PED_NOISE << endl;


// load efficiency correction factors
Expand Down Expand Up @@ -68,11 +80,6 @@ fdc_config_t::fdc_config_t(const std::shared_ptr<const JEvent>& event)
channel_efficiencies.push_back( new_strip_efficiencies[2*chamber] );
}
}

FDC_EFFVSDOCA_PAR[0] = 0.999;
FDC_EFFVSDOCA_PAR[1] = 3.75e-4;
FDC_EFFVSDOCA_PAR[2] = 0.506;
FDC_EFFVSDOCA_PAR[3] = 3.75e-2;
}


Expand All @@ -83,10 +90,53 @@ void FDCSmearer::SmearEvent(hddm_s::HDDM *record)
{
double t_max = config->TRIGGER_LOOKBACK_TIME + fdc_config->FDC_TIME_WINDOW;
double threshold = fdc_config->FDC_THRESHOLD_FACTOR * fdc_config->FDC_PED_NOISE; // for sparsification
if (fdc_config->FDC_STRIP_THR_SET >= 0) threshold = fdc_config->FDC_STRIP_THR_SET;
double wire_thr_kev = (fdc_config->FDC_WIRE_THR_KEV_SET >= 0) ? fdc_config->FDC_WIRE_THR_KEV_SET : fdc_config->FDC_THRESH_KEV;
double gain_global = fdc_config->FDC_GAIN_SCALE; // gas gain relative to the simulation: scales the avalanche charge seen on wires and strips

hddm_s::FdcChamberList chambers = record->getFdcChambers();
hddm_s::FdcChamberList::iterator iter;
for (iter = chambers.begin(); iter != chambers.end(); ++iter) {
int pkg = (iter->getModule() - 1) / 2; // FDC package 0-3 (3 planes per module, 2 modules per package)
double gain = gain_global * ((pkg >= 0 && pkg < 4) ? fdc_config->FDC_GAIN_PKG[pkg] : 1.);

// drift-time fluctuation common to the wire hit and the strip hits induced by the same track in this chamber
// (the strips see the charge induced by the wire avalanche, so the drift delay is the same for both); one shift per itrack
std::map<int, double> drift_shift_ns;
auto drift_shift = [&](int itrack) {
if (!config->SMEAR_HITS || fdc_config->FDC_DRIFT_SMEAR_NS <= 0) return 0.;
auto it = drift_shift_ns.find(itrack);
if (it == drift_shift_ns.end())
it = drift_shift_ns.emplace(itrack, gDRandom.SampleGaussian(fdc_config->FDC_DRIFT_SMEAR_NS)).first;
return it->second;
};

// signal attenuation vs the track's doca to its wire in this chamber (charge lost at the cell edge), per itrack, from the anode truth hits
std::map<int, double> atten;
auto atten_of_doca = [&](double d) { return std::max(0.05, 1. - fdc_config->FDC_DOCA_ATTEN * d * d); };
if (fdc_config->FDC_DOCA_ATTEN > 0) {
hddm_s::FdcAnodeWireList awires = iter->getFdcAnodeWires();
for (auto aw = awires.begin(); aw != awires.end(); ++aw) {
hddm_s::FdcAnodeTruthHitList ath = aw->getFdcAnodeTruthHits();
for (auto at = ath.begin(); at != ath.end(); ++at) {
double a = atten_of_doca(at->getD());
auto it = atten.find(at->getItrack());
if (it == atten.end() || a > it->second) atten[at->getItrack()] = a; // the wire with the closest approach
}
}
}
auto atten_of_track = [&](int itrack) { auto it = atten.find(itrack); return (it == atten.end()) ? 1. : it->second; };
// avalanche gain fluctuation, one factor per track in this chamber (the strips see the charge induced by the same wire avalanche)
std::map<int, double> fluct;
auto fluct_of_track = [&](int itrack) {
if (fdc_config->FDC_GAIN_SIGMA <= 0) return 1.;
auto it = fluct.find(itrack);
if (it == fluct.end()) {
double s = fdc_config->FDC_GAIN_SIGMA;
it = fluct.emplace(itrack, exp(gDRandom.SampleGaussian(s) - 0.5 * s * s)).first;
}
return it->second;
};

// Add pedestal noise to strip charge data
hddm_s::FdcCathodeStripList strips = iter->getFdcCathodeStrips();
Expand All @@ -103,8 +153,8 @@ void FDCSmearer::SmearEvent(hddm_s::HDDM *record)
&& !gDRandom.DecideToAcceptHit(fdc_config->GetEfficiencyCorrectionFactor(siter)))
continue;

double q = titer->getQ();
double t = titer->getT();
double q = titer->getQ() * gain * atten_of_track(titer->getItrack()) * fluct_of_track(titer->getItrack());
double t = titer->getT() + drift_shift(titer->getItrack());
if(config->SMEAR_HITS) {
q += gDRandom.SampleGaussian(fdc_config->FDC_PED_NOISE);
t += gDRandom.SampleGaussian(fdc_config->FDC_TDRIFT_SIGMA)*1.0e9;
Expand All @@ -129,16 +179,17 @@ void FDCSmearer::SmearEvent(hddm_s::HDDM *record)
hddm_s::FdcAnodeTruthHitList thits = witer->getFdcAnodeTruthHits();
hddm_s::FdcAnodeTruthHitList::iterator titer;
for (titer = thits.begin(); titer != thits.end(); ++titer) {
// anode discriminator threshold, applied to the energy deposited on the wire
// (dE in GeV, summed by hdgeant4 over the track segments it combined into this hit)
double att = (fdc_config->FDC_DOCA_ATTEN > 0) ? atten_of_doca(titer->getD()) : 1.;
if (titer->getDE() * gain * att * fluct_of_track(titer->getItrack()) * 1e6 < wire_thr_kev)
continue;
// correct simulation efficiencies
if (config->APPLY_EFFICIENCY_CORRECTIONS
&& !gDRandom.DecideToAcceptHit(fdc_config->GetEfficiencyCorrectionFactor(witer)))
continue;
double doca = titer->getD();
if (config->APPLY_EFFICIENCY_CORRECTIONS
&& !gDRandom.DecideToAcceptHit(fdc_config->GetEfficiencyVsDOCA(doca)))
continue;

double t = titer->getT();
double t = titer->getT() + drift_shift(titer->getItrack());
if(config->SMEAR_HITS) {
t += gDRandom.SampleGaussian(fdc_config->FDC_TDRIFT_SIGMA)*1.0e9;
}
Expand Down
23 changes: 8 additions & 15 deletions src/programs/Simulation/mcsmear/FDCSmearer.h
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,14 @@ class fdc_config_t
double FDC_TIME_WINDOW;
double FDC_THRESH_KEV; // fdc anode discriminator threshold

static bool FDC_EFFVSDOCA;
double FDC_EFFVSDOCA_PAR[4];
static double FDC_DRIFT_SMEAR_NS; // FDC:DRIFT_SMEAR_NS, Gaussian drift-time smear (ns) common to the wire and strip hits of one track in one chamber (0 = off)
static double FDC_CATHODE_SIGMA_SET; // FDC:CATHODE_SIGMA, override of FDC_CATHODE_SIGMA in microns (< 0 = use CCDB)
static double FDC_GAIN_SCALE; // FDC:GAIN_SCALE, gas-gain factor applied to the strip charge q and the wire dE before the thresholds (1 = off)
static double FDC_GAIN_PKG[4]; // FDC:GAIN_PKG1..4, extra gain factor per FDC package (modules 1-2, 3-4, 5-6, 7-8), times FDC_GAIN_SCALE (1 = off)
static double FDC_STRIP_THR_SET; // FDC:STRIP_THRESHOLD, strip charge threshold in pC (< 0 = FDC_THRESHOLD_FACTOR * FDC_PED_NOISE from CCDB)
static double FDC_GAIN_SIGMA; // FDC:GAIN_SIGMA s, avalanche gain fluctuation: one log-normal factor exp(s G - s^2/2) per chamber and track, shared by its wire hit and the strips it induces (0 = off)
static double FDC_DOCA_ATTEN; // FDC:DOCA_ATTEN a, signal (wire dE and the track's strip q in the chamber) scaled by max(0.05, 1 - a doca^2), doca in cm (0 = off)
static double FDC_WIRE_THR_KEV_SET; // FDC:WIRE_THRESHOLD_KEV, anode threshold on the gain-scaled wire dE in keV (< 0 = FDC_THRESH_KEV from CCDB)

vector< vector<double > > channel_efficiencies;

Expand Down Expand Up @@ -49,19 +55,6 @@ class fdc_config_t

return channel_efficiencies.at(gPlane).at(element-1);
}

double GetEfficiencyVsDOCA(double doca) {
// Introduce a DOCA-dependent efficiency that cuts off the hit efficiency
// for tracks that cross through the cell at a large distance from the wire.
double eff = FDC_EFFVSDOCA_PAR[0] - FDC_EFFVSDOCA_PAR[1] *
(1 / (pow(FDC_EFFVSDOCA_PAR[2] - doca, 2) +
pow(FDC_EFFVSDOCA_PAR[3], 2))
- 2 * FDC_EFFVSDOCA_PAR[2] * doca /
pow(pow(FDC_EFFVSDOCA_PAR[2], 2) +
pow(FDC_EFFVSDOCA_PAR[3], 2), 2)
);
return (!FDC_EFFVSDOCA)? 1 : (eff < 0)? 0 : eff*eff;
}
};


Expand Down
24 changes: 17 additions & 7 deletions src/programs/Simulation/mcsmear/MyProcessor.cc
Original file line number Diff line number Diff line change
Expand Up @@ -159,13 +159,23 @@ void MyProcessor::Init(void)
bzero(&input_file_mutex_last_owner, sizeof(pthread_t));
bzero(&smearer_mutex_last_owner, sizeof(pthread_t));

// By default, the empirical fdc DOCA-dependent efficiency
// is applied to fdc wire hits in FDCSmearer. This can be
// disabled by this commandline option.
app->SetDefaultParameter("FDC:USE_EFFVSDOCA", fdc_config_t::FDC_EFFVSDOCA,
"Turn on/off fdc doca-dependent wire hit efficiency,"
" \"0\"=no efficiency correction,"
" \"1\"=standard efficiency correction (default)");
app->SetDefaultParameter("FDC:DRIFT_SMEAR_NS", fdc_config_t::FDC_DRIFT_SMEAR_NS,
"Gaussian drift-time smear in ns, common to the wire and strip hits of one track in one chamber; 0 = off");
app->SetDefaultParameter("FDC:CATHODE_SIGMA", fdc_config_t::FDC_CATHODE_SIGMA_SET,
"FDC_CATHODE_SIGMA override in microns, sets the strip pedestal noise; < 0 = CCDB");
app->SetDefaultParameter("FDC:GAIN_SCALE", fdc_config_t::FDC_GAIN_SCALE,
"FDC gas-gain factor on the strip charge and wire dE before the thresholds; 1 = off");
for (int i = 0; i < 4; i++)
app->SetDefaultParameter("FDC:GAIN_PKG" + std::to_string(i + 1), fdc_config_t::FDC_GAIN_PKG[i],
"FDC extra gain factor for package " + std::to_string(i + 1) + " (times FDC:GAIN_SCALE); 1 = off");
app->SetDefaultParameter("FDC:GAIN_SIGMA", fdc_config_t::FDC_GAIN_SIGMA,
"FDC avalanche gain fluctuation (log-normal sigma), one factor per chamber and track for wire and strips; 0 = off");
app->SetDefaultParameter("FDC:DOCA_ATTEN", fdc_config_t::FDC_DOCA_ATTEN,
"FDC signal attenuation a: wire dE and strip q scaled by max(0.05, 1 - a doca^2) (doca in cm); 0 = off");
app->SetDefaultParameter("FDC:STRIP_THRESHOLD", fdc_config_t::FDC_STRIP_THR_SET,
"FDC strip charge threshold in pC; < 0 = FDC_THRESHOLD_FACTOR * FDC_PED_NOISE (CCDB)");
app->SetDefaultParameter("FDC:WIRE_THRESHOLD_KEV", fdc_config_t::FDC_WIRE_THR_KEV_SET,
"FDC anode threshold on the gain-scaled wire dE in keV; < 0 = FDC_THRESH_KEV (CCDB)");

return; //NOERROR;
}
Expand Down
Loading