diff --git a/src/programs/Simulation/mcsmear/FDCSmearer.cc b/src/programs/Simulation/mcsmear/FDCSmearer.cc index 96b1f83cd..2dbf8041b 100644 --- a/src/programs/Simulation/mcsmear/FDCSmearer.cc +++ b/src/programs/Simulation/mcsmear/FDCSmearer.cc @@ -7,7 +7,14 @@ #include -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) @@ -40,7 +47,12 @@ fdc_config_t::fdc_config_t(const std::shared_ptr& 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 @@ -68,11 +80,6 @@ fdc_config_t::fdc_config_t(const std::shared_ptr& 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; } @@ -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 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 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 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(); @@ -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; @@ -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; } diff --git a/src/programs/Simulation/mcsmear/FDCSmearer.h b/src/programs/Simulation/mcsmear/FDCSmearer.h index d710cb636..24858421e 100644 --- a/src/programs/Simulation/mcsmear/FDCSmearer.h +++ b/src/programs/Simulation/mcsmear/FDCSmearer.h @@ -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 > channel_efficiencies; @@ -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; - } }; diff --git a/src/programs/Simulation/mcsmear/MyProcessor.cc b/src/programs/Simulation/mcsmear/MyProcessor.cc index b05db5af6..1ff0ee1f2 100644 --- a/src/programs/Simulation/mcsmear/MyProcessor.cc +++ b/src/programs/Simulation/mcsmear/MyProcessor.cc @@ -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; }