From 815e385c03d985fa55ea34db6229f10b689ea2e0 Mon Sep 17 00:00:00 2001 From: Richard Jones Date: Mon, 2 Nov 2026 21:50:34 -0500 Subject: [PATCH 1/2] mcsmear: write the FDC strip peak amplitude (fdcDigihit peakAmp) and carry it through the random-trigger merge [rtj] In real data the FDC cathode strips are read out by fADC125s, and the reconstruction uses the pulse peak amplitude (DFDCHit::pulse_height, with q = a_scale * pulse_height). mcsmear wrote only the integrated charge q, so for MC the reconstruction fell back on q where the data use the peak. FDCSmearer now adds an fdcDigihit to every smeared strip hit, with peakAmp = q / FDC_ADC_ASCALE (/FDC/digi_scales), i.e. on the data scale. If FDC_ADC_ASCALE cannot be read from CCDB no fdcDigihit is written and the behaviour is as before. The FDC strip merge with random triggers (hddm_s_merger, the earliest-hit rule of PR #420) now applies the same rule to the peak amplitude as to q: the earliest random hit replaces the MC hit (time, q and peakAmp) if it is earlier, both are summed if they fall in the same 8 ns sample, and a later random hit is ignored. A hit without an fdcDigihit contributes q, which is what the reconstruction would use for it. The merge also returns immediately if the random strip has no hits, instead of dereferencing src.begin() of an empty list. Validated in the GlueX rho0 data/MC study of 2026-10/11: with this mcsmear (and the hdgeant4 FDC fixes of rjones30/HDGeant4 PR #37) the change of the rho0 yield between recon-2018_08 and recon-2019_11 on the same data is reproduced by the MC to 0.3% (70 chunks of 10k events, runs 51048-51052). Co-Authored-By: Claude Opus 5.5 (cherry picked from commit c362337e1b8b6ef61fc8f6447b6c6405f1252489) --- src/programs/Simulation/mcsmear/FDCSmearer.cc | 14 ++++++++++ src/programs/Simulation/mcsmear/FDCSmearer.h | 1 + .../Simulation/mcsmear/hddm_s_merger.cc | 27 +++++++++++++++++++ 3 files changed, 42 insertions(+) diff --git a/src/programs/Simulation/mcsmear/FDCSmearer.cc b/src/programs/Simulation/mcsmear/FDCSmearer.cc index 96b1f83cd..c99ddc13d 100644 --- a/src/programs/Simulation/mcsmear/FDCSmearer.cc +++ b/src/programs/Simulation/mcsmear/FDCSmearer.cc @@ -37,6 +37,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 @@ -113,6 +122,11 @@ void FDCSmearer::SmearEvent(hddm_s::HDDM *record) 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..d9816ba59 100644 --- a/src/programs/Simulation/mcsmear/FDCSmearer.h +++ b/src/programs/Simulation/mcsmear/FDCSmearer.h @@ -18,6 +18,7 @@ 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; double FDC_EFFVSDOCA_PAR[4]; 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. From 14a1e5f079b558ab40dac88efdec9ac4a1a1eafb Mon Sep 17 00:00:00 2001 From: Richard Jones Date: Mon, 9 Nov 2026 14:58:26 -0500 Subject: [PATCH 2/2] mcsmear: FDC cathode strip readout thresholds on the raw amplitude q/a_scale, per strip from CCDB [rtj] A strip hit is dropped if raw_amp = q / FDC_ADC_ASCALE is below its readout threshold in ADC counts. The thresholds are read per strip from the CCDB tables /FDC/packageN/strip_thresholds (12 rows per package, in the row order of strip_mc_efficiency); where a table is absent every strip of that package uses 50 counts. FDC:RAW_AMP_THRESHOLD >= 0 overrides the tables with one value for every strip; 0 is no cut, since strips with q <= 0 are dropped by the sparsification threshold anyway. The cut draws no random numbers. No cut is applied if FDC/digi_scales FDC_ADC_ASCALE is missing. Requested by N. Jarvis on PR #422 as a separate PR; needs the peakAmp PR (FDC_ADC_ASCALE). Supersedes bd0cb53, whose single value defaulted to off. Co-Authored-By: Claude Opus 5.5 --- src/programs/Simulation/mcsmear/FDCSmearer.cc | 30 +++++++++++++++++++ src/programs/Simulation/mcsmear/FDCSmearer.h | 15 ++++++++++ .../Simulation/mcsmear/MyProcessor.cc | 6 ++++ 3 files changed, 51 insertions(+) diff --git a/src/programs/Simulation/mcsmear/FDCSmearer.cc b/src/programs/Simulation/mcsmear/FDCSmearer.cc index c99ddc13d..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) @@ -78,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; @@ -118,6 +144,10 @@ 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); diff --git a/src/programs/Simulation/mcsmear/FDCSmearer.h b/src/programs/Simulation/mcsmear/FDCSmearer.h index d9816ba59..cf1993931 100644 --- a/src/programs/Simulation/mcsmear/FDCSmearer.h +++ b/src/programs/Simulation/mcsmear/FDCSmearer.h @@ -21,9 +21,13 @@ class fdc_config_t 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 @@ -38,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,"