diff --git a/src/programs/Simulation/mcsmear/FDCSmearer.cc b/src/programs/Simulation/mcsmear/FDCSmearer.cc index 96b1f83cd..cb0af0c0b 100644 --- a/src/programs/Simulation/mcsmear/FDCSmearer.cc +++ b/src/programs/Simulation/mcsmear/FDCSmearer.cc @@ -8,6 +8,8 @@ #include bool fdc_config_t::FDC_EFFVSDOCA(1); +double fdc_config_t::FDC_RAW_AMP_THR(-1.); +const double fdc_config_t::FDC_RAW_AMP_THR_DEFAULT(50.); //----------- // fdc_config_t (constructor) @@ -37,6 +39,15 @@ fdc_config_t::fdc_config_t(const std::shared_ptr& event) FDC_THRESH_KEV = fdcparms["FDC_THRESH_KEV"]; } + // scale for the strip peak amplitude written in fdcDigihit (0 = do not write fdcDigihit) + FDC_ADC_ASCALE = 0.; + map digi_scales; + if(DEvent::GetCalib(event, "FDC/digi_scales", digi_scales) || digi_scales.find("FDC_ADC_ASCALE") == digi_scales.end()) + jerr << "Problem loading FDC/digi_scales FDC_ADC_ASCALE from CCDB! fdcDigihit peakAmp will not be written." << endl; + else + FDC_ADC_ASCALE = digi_scales["FDC_ADC_ASCALE"]; + jout << "FDC strip peakAmp = q / FDC_ADC_ASCALE, FDC_ADC_ASCALE = " << FDC_ADC_ASCALE << endl; + // 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 @@ -69,6 +80,30 @@ fdc_config_t::fdc_config_t(const std::shared_ptr& event) } } + // readout thresholds of the cathode strips on raw_amp = q / FDC_ADC_ASCALE, in ADC counts, per strip from + // /FDC/packageN/strip_thresholds (12 rows, as strip_mc_efficiency); FDC_RAW_AMP_THR_DEFAULT where a table is absent + for(int package=1; package<=4; package++) { + vector< vector > new_thresholds; + char ccdb_str[100]; + sprintf(ccdb_str, "/FDC/package%d/strip_thresholds", package); + bool have_table = false; + if (FDC_RAW_AMP_THR < 0) { + have_table = !DEvent::GetCalib(event, ccdb_str, new_thresholds) && new_thresholds.size() == 12; + if (have_table) + jout << "FDC strip raw-amplitude thresholds of package " << package << " from " << ccdb_str << endl; + else + jout << "FDC strip raw-amplitude thresholds of package " << package << ": no " << ccdb_str + << ", default " << FDC_RAW_AMP_THR_DEFAULT << " ADC counts" << endl; + } + for(int chamber=0; chamber<6; chamber++) { + channel_raw_amp_thresholds.push_back( have_table? new_thresholds[2*chamber+1] : vector() ); + channel_raw_amp_thresholds.push_back( vector() ); + channel_raw_amp_thresholds.push_back( have_table? new_thresholds[2*chamber] : vector() ); + } + } + if (FDC_RAW_AMP_THR >= 0) + jout << "FDC strip raw-amplitude threshold FDC:RAW_AMP_THRESHOLD = " << FDC_RAW_AMP_THR << " ADC counts for every strip" << endl; + FDC_EFFVSDOCA_PAR[0] = 0.999; FDC_EFFVSDOCA_PAR[1] = 3.75e-4; FDC_EFFVSDOCA_PAR[2] = 0.506; @@ -109,10 +144,19 @@ void FDCSmearer::SmearEvent(hddm_s::HDDM *record) q += gDRandom.SampleGaussian(fdc_config->FDC_PED_NOISE); t += gDRandom.SampleGaussian(fdc_config->FDC_TDRIFT_SIGMA)*1.0e9; } + // readout threshold on the raw amplitude raw_amp = q / FDC_ADC_ASCALE, in ADC counts + if (fdc_config->FDC_ADC_ASCALE > 0 + && q / fdc_config->FDC_ADC_ASCALE < fdc_config->GetRawAmpThreshold(siter)) + continue; if (q > threshold && t > config->TRIGGER_LOOKBACK_TIME && t < t_max) { hddm_s::FdcCathodeHitList hits = siter->addFdcCathodeHits(); hits().setQ(q); hits().setT(t); + // strip peak amplitude on the data scale, read by DEventSourceHDDM as DFDCHit::pulse_height + if (fdc_config->FDC_ADC_ASCALE > 0) { + hddm_s::FdcDigihitList digihit = hits().addFdcDigihits(); + digihit().setPeakAmp(q / fdc_config->FDC_ADC_ASCALE); + } } } diff --git a/src/programs/Simulation/mcsmear/FDCSmearer.h b/src/programs/Simulation/mcsmear/FDCSmearer.h index d710cb636..cf1993931 100644 --- a/src/programs/Simulation/mcsmear/FDCSmearer.h +++ b/src/programs/Simulation/mcsmear/FDCSmearer.h @@ -18,11 +18,16 @@ class fdc_config_t //double FDC_HIT_DROP_FRACTION = 0.0; // 1000.0E-9 double FDC_TIME_WINDOW; double FDC_THRESH_KEV; // fdc anode discriminator threshold + double FDC_ADC_ASCALE; // /FDC/digi_scales FDC_ADC_ASCALE; strip peakAmp = q / FDC_ADC_ASCALE, as in data (DFDCHit_factory: q = a_scale * pulse_height) static bool FDC_EFFVSDOCA; + static double FDC_RAW_AMP_THR; // FDC:RAW_AMP_THRESHOLD in ADC counts; < 0 (default): per strip from CCDB /FDC/packageN/strip_thresholds, + // FDC_RAW_AMP_THR_DEFAULT where a table is absent; >= 0: this value for every strip (0: no cut) + static const double FDC_RAW_AMP_THR_DEFAULT; double FDC_EFFVSDOCA_PAR[4]; vector< vector > channel_efficiencies; + vector< vector > channel_raw_amp_thresholds; // same row order as channel_efficiencies; empty rows for the anode wires double GetEfficiencyCorrectionFactor(hddm_s::FdcCathodeStripList::iterator &siter) { // cathode strips @@ -37,6 +42,17 @@ class fdc_config_t return channel_efficiencies.at(gPlane).at(element-1); } + double GetRawAmpThreshold(hddm_s::FdcCathodeStripList::iterator &siter) { + // readout threshold of a cathode strip on raw_amp = q / FDC_ADC_ASCALE, in ADC counts + if (FDC_RAW_AMP_THR >= 0) + return FDC_RAW_AMP_THR; + int gPlane = 9*(siter->getModule()-1) + 3*(siter->getLayer()-1) + + (siter->getPlane()-1); // starts counting at 0 + int element = siter->getStrip(); + const vector &row = channel_raw_amp_thresholds.at(gPlane); + return (element >= 1 && element <= (int)row.size())? row[element-1] : FDC_RAW_AMP_THR_DEFAULT; + } + double GetEfficiencyCorrectionFactor(hddm_s::FdcAnodeWireList::iterator &witer) { // anode wires int gPlane = 9*(witer->getModule()-1) + 3*(witer->getLayer()-1) + 1; // starts counting at 0 diff --git a/src/programs/Simulation/mcsmear/MyProcessor.cc b/src/programs/Simulation/mcsmear/MyProcessor.cc index b05db5af6..bb5c782ea 100644 --- a/src/programs/Simulation/mcsmear/MyProcessor.cc +++ b/src/programs/Simulation/mcsmear/MyProcessor.cc @@ -162,6 +162,12 @@ void MyProcessor::Init(void) // By default, the empirical fdc DOCA-dependent efficiency // is applied to fdc wire hits in FDCSmearer. This can be // disabled by this commandline option. + // FDC cathode strip readout threshold on raw_amp = q / FDC_ADC_ASCALE (ADC counts). By default the + // per-strip values of /FDC/packageN/strip_thresholds are used, or 50 where a table is absent. + app->SetDefaultParameter("FDC:RAW_AMP_THRESHOLD", fdc_config_t::FDC_RAW_AMP_THR, + "FDC strip dropped if raw_amp = q / FDC_ADC_ASCALE is below this, in ADC counts;" + " < 0 = per strip from CCDB /FDC/packageN/strip_thresholds, 50 where absent (default);" + " >= 0 = this value for every strip (0 = no cut)"); app->SetDefaultParameter("FDC:USE_EFFVSDOCA", fdc_config_t::FDC_EFFVSDOCA, "Turn on/off fdc doca-dependent wire hit efficiency," " \"0\"=no efficiency correction," diff --git a/src/programs/Simulation/mcsmear/hddm_s_merger.cc b/src/programs/Simulation/mcsmear/hddm_s_merger.cc index 69be7ff16..91a388510 100644 --- a/src/programs/Simulation/mcsmear/hddm_s_merger.cc +++ b/src/programs/Simulation/mcsmear/hddm_s_merger.cc @@ -960,9 +960,24 @@ hddm_s::FdcCathodeStripList &operator+=(hddm_s::FdcCathodeStripList &dst, return dst; } +// Peak amplitude of an FDC strip hit, as the reconstruction reads it: the +// fdcDigihit peakAmp if there is one, otherwise q. If the hit has no +// fdcDigihit, one holding q is added so that the merged amplitude can be +// stored in it. +static double fdc_strip_peakamp(hddm_s::FdcCathodeHit &hit) +{ + if (hit.getFdcDigihits().size() > 0) + return hit.getFdcDigihit().getPeakAmp(); + hddm_s::FdcDigihitList digihit = hit.addFdcDigihits(); + digihit().setPeakAmp(hit.getQ()); + return hit.getQ(); +} + hddm_s::FdcCathodeHitList &operator+=(hddm_s::FdcCathodeHitList &dst, hddm_s::FdcCathodeHitList &src) { + if (src.size() == 0) + return dst; // There should be no more than one hit in dst - assume this is the case. (otherwise, would have to pick the earliest) @@ -970,6 +985,8 @@ hddm_s::FdcCathodeHitList &operator+=(hddm_s::FdcCathodeHitList &dst, // If it is before the original hit, replace the original hit's time and peak height. // If it is at the same time sample as the original hit, sum the pulse height. // If it is after the original hit, leave the original hit's time and peak height alone. + // The strip peak amplitude (fdcDigihit peakAmp, read by the reconstruction as the + // pulse height) follows the same rule as q. hddm_s::FdcCathodeHitList::iterator iter; @@ -990,11 +1007,17 @@ hddm_s::FdcCathodeHitList &operator+=(hddm_s::FdcCathodeHitList &dst, if (t > fdc_strips_readout_window_ns) return dst; // random is too late double newQ = iter_firsthit->getQ(); + bool newA_exists = (iter_firsthit->getFdcDigihits().size() > 0); + double newA = (newA_exists)? iter_firsthit->getFdcDigihit().getPeakAmp() : newQ; if (dst.size() == 0) { // no hits exist, add the random dst.add(1, -1); dst(0).setT(t); dst(0).setQ(newQ); + if (newA_exists) { + hddm_s::FdcDigihitList digihit = dst(0).addFdcDigihits(); + digihit().setPeakAmp(newA); + } } else { @@ -1009,13 +1032,17 @@ hddm_s::FdcCathodeHitList &operator+=(hddm_s::FdcCathodeHitList &dst, if (src_hitsample == dst_hitsample) { // same sample, add pulse heights + double origA = fdc_strip_peakamp(dst(0)); if (t < origT) dst(0).setT(t); dst(0).setQ(origQ + newQ); + dst(0).getFdcDigihit().setPeakAmp(origA + newA); } else if (src_hitsample < dst_hitsample) { // random arrives earlier, replace hit time and pulse height + fdc_strip_peakamp(dst(0)); dst(0).setT(t); dst(0).setQ(newQ); + dst(0).getFdcDigihit().setPeakAmp(newA); } // if the random hit is after the original one, it is too late, ignore it.