From 815e385c03d985fa55ea34db6229f10b689ea2e0 Mon Sep 17 00:00:00 2001 From: Richard Jones Date: Mon, 2 Nov 2026 21:50:34 -0500 Subject: [PATCH] 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.