From 57010d15f8512b9eb7708493752632c71af71e1f Mon Sep 17 00:00:00 2001 From: "poppi@baltig.infn.it" Date: Tue, 30 Jun 2026 18:40:33 +0200 Subject: [PATCH 1/5] Addition of some WireMod diagnostic tools, addition of some flags (usage of SignalROIF, rounding of floats, fixed scaling factors to 1 at al. ) --- sbncode/HitFinder/ChannelROIToWire_module.cc | 29 ++- sbncode/WireMod/Utility/CMakeLists.txt | 1 + sbncode/WireMod/Utility/WireModUtility.cc | 204 ++++++++++++++++--- sbncode/WireMod/Utility/WireModUtility.hh | 10 +- 4 files changed, 202 insertions(+), 42 deletions(-) diff --git a/sbncode/HitFinder/ChannelROIToWire_module.cc b/sbncode/HitFinder/ChannelROIToWire_module.cc index 42d304537..497ffedcb 100644 --- a/sbncode/HitFinder/ChannelROIToWire_module.cc +++ b/sbncode/HitFinder/ChannelROIToWire_module.cc @@ -51,6 +51,7 @@ class ChannelROIToWire : public art::EDProducer std::vector fWireModuleLabelVec; ///< vector of modules that made digits std::vector fOutInstanceLabelVec; ///< The output instance labels to apply bool fDiagnosticOutput; ///< secret diagnostics flag + bool fUseFloatPrecision; ///< use SignalROIF() to preserve sub-ADC precision (no rounding) size_t fEventCount; ///< count of event processed const geo::WireReadoutGeom* fChannelMapAlg = &art::ServiceHandle()->Get(); @@ -77,6 +78,7 @@ void ChannelROIToWire::reconfigure(fhicl::ParameterSet const& pset) fWireModuleLabelVec = pset.get>("WireModuleLabelVec", std::vector()={"decon1droi"}); fOutInstanceLabelVec = pset.get> ("OutInstanceLabelVec", {"PHYSCRATEDATA"}); fDiagnosticOutput = pset.get< bool >("DiagnosticOutput", false); + fUseFloatPrecision = pset.get< bool >("UseFloatPrecision", false); if (fWireModuleLabelVec.size() != fOutInstanceLabelVec.size()) { @@ -131,19 +133,24 @@ void ChannelROIToWire::produce(art::Event& evt) // Create an ROI vector for output recob::Wire::RegionsOfInterest_t ROIVec; - - // Loop through the ROIs for this channel - const recob::ChannelROI::RegionsOfInterest_t& channelROIs = channelROI.SignalROI(); - for(const auto& range : channelROIs.get_ranges()) + if (fUseFloatPrecision) { - size_t startTick = range.begin_index(); - - std::vector dataVec(range.data().size()); - - for(size_t binIdx = 0; binIdx < range.data().size(); binIdx++) dataVec[binIdx] = std::round(range.data()[binIdx] / ADCScaleFactor); - - ROIVec.add_range(startTick, std::move(dataVec)); + const recob::ChannelROI::RegionsOfInterest_f& channelROIsF = channelROI.SignalROIF(); + for(const auto& range : channelROIsF.get_ranges()) + ROIVec.add_range(range.begin_index(), range.data()); + } + else + { + const recob::ChannelROI::RegionsOfInterest_t& channelROIs = channelROI.SignalROI(); + for(const auto& range : channelROIs.get_ranges()) + { + size_t startTick = range.begin_index(); + std::vector dataVec(range.data().size()); + for(size_t binIdx = 0; binIdx < range.data().size(); binIdx++) + dataVec[binIdx] = std::round(range.data()[binIdx] / ADCScaleFactor); + ROIVec.add_range(startTick, std::move(dataVec)); + } } wireCol->push_back(recob::WireCreator(std::move(ROIVec),channel,view).move()); diff --git a/sbncode/WireMod/Utility/CMakeLists.txt b/sbncode/WireMod/Utility/CMakeLists.txt index 8955de261..3c441620f 100644 --- a/sbncode/WireMod/Utility/CMakeLists.txt +++ b/sbncode/WireMod/Utility/CMakeLists.txt @@ -17,6 +17,7 @@ cet_make_library( lardataobj::AnalysisBase lardataobj::RawData lardataobj::RecoBase + sbnobj::ICARUS_TPC lardataalg::DetectorInfo lardata::RecoObjects lardata::Utilities diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 8ca184dc5..3b583f048 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -43,6 +43,46 @@ sys::WireModUtility::ROIProperties_t sys::WireModUtility::CalcROIProperties(reco return roi_vals; } +sys::WireModUtility::ROIProperties_t sys::WireModUtility::CalcROIProperties(recob::ChannelROI const& chanROI, size_t const& roi_idx) +{ + // SignalROIF() returns by value, so use SignalROI() and convert short->float manually. + auto const& roi_short = chanROI.SignalROI().get_ranges()[roi_idx]; + const float scale = 1.0f / static_cast(chanROI.ADCScaleFactor()); + + ROIProperties_t roi_vals; + roi_vals.channel = chanROI.Channel(); + roi_vals.view = wireReadout->View(chanROI.Channel()); + roi_vals.begin = roi_short.begin_index(); + roi_vals.end = roi_short.end_index(); + roi_vals.center = 0; + roi_vals.total_q = 0; + roi_vals.sigma = 0; + + auto const& short_data = roi_short.data(); + for (size_t i_t = 0; i_t < short_data.size(); ++i_t) + { + float val = static_cast(short_data[i_t]) * scale; + roi_vals.center += val * (i_t + roi_vals.begin); + roi_vals.total_q += val; + } + roi_vals.center = roi_vals.center / roi_vals.total_q; + + for (size_t i_t = 0; i_t < short_data.size(); ++i_t) + { + float val = static_cast(short_data[i_t]) * scale; + roi_vals.sigma += val * (i_t + roi_vals.begin - roi_vals.center) + * (i_t + roi_vals.begin - roi_vals.center); + } + roi_vals.sigma = std::sqrt(roi_vals.sigma / roi_vals.total_q); + if (roi_vals.end - roi_vals.begin == 1) + { + roi_vals.center += 0.5; + roi_vals.sigma = 0.5; + } + + return roi_vals; +} + //--- GetTargetROIs --- std::vector> sys::WireModUtility::GetTargetROIs(sim::SimEnergyDeposit const& shifted_edep, double offset) { @@ -204,6 +244,78 @@ void sys::WireModUtility::FillROIMatchedHitMap(std::vector const& hi } } +//--- FillROIMatchedEdepMap (ChannelROI overload) --- +void sys::WireModUtility::FillROIMatchedEdepMap(std::vector const& edepVec, std::vector const& chanROIVec, double offset) +{ + ROIMatchedEdepMap.clear(); + + std::unordered_map chanROIChannelMap; + for (size_t i_c = 0; i_c < chanROIVec.size(); ++i_c) + chanROIChannelMap[chanROIVec[i_c].Channel()] = i_c; + + for (size_t i_e = 0; i_e < edepVec.size(); ++i_e) + { + std::vector> target_rois = GetTargetROIs(edepVec[i_e], offset); + + for (auto const& target_roi : target_rois) + { + if (chanROIChannelMap.find(target_roi.first) == chanROIChannelMap.end()) + continue; + + auto const& target_chanROI = chanROIVec.at(chanROIChannelMap[target_roi.first]); + + // SignalROIF() returns by value, so use SignalROI() instead. + auto const& roiShort = target_chanROI.SignalROI(); + if (not roiShort.is_valid() || + roiShort.empty() || + roiShort.n_ranges() == 0 || + roiShort.size() <= target_roi.second || + roiShort.is_void(target_roi.second) ) + continue; + + auto range_number = roiShort.find_range_iterator(target_roi.second) - roiShort.begin_range(); + + ROIMatchedEdepMap[std::make_pair(target_chanROI.Channel(), range_number)].push_back(i_e); + } + } +} + +//--- FillROIMatchedHitMap (ChannelROI overload) --- +void sys::WireModUtility::FillROIMatchedHitMap(std::vector const& hitVec, std::vector const& chanROIVec) +{ + ROIMatchedHitMap.clear(); + + std::unordered_map chanROIChannelMap; + for (size_t i_c = 0; i_c < chanROIVec.size(); ++i_c) + chanROIChannelMap[chanROIVec[i_c].Channel()] = i_c; + + for (size_t i_h = 0; i_h < hitVec.size(); ++i_h) + { + std::vector> target_rois = GetHitTargetROIs(hitVec[i_h]); + + for (auto const& target_roi : target_rois) + { + if (chanROIChannelMap.find(target_roi.first) == chanROIChannelMap.end()) + continue; + + auto const& target_chanROI = chanROIVec.at(chanROIChannelMap[target_roi.first]); + + // SignalROIF() returns by value, so use SignalROI() instead. + auto const& roiShort = target_chanROI.SignalROI(); + if (not roiShort.is_valid() || + roiShort.empty() || + roiShort.n_ranges() == 0 || + roiShort.size() <= target_roi.second || + roiShort.is_void(target_roi.second) ) + continue; + + auto range_number = roiShort.find_range_iterator(target_roi.second) - roiShort.begin_range(); + + ROIMatchedHitMap[std::make_pair(target_chanROI.Channel(), range_number)].push_back(i_h); + } + } +} + //--- CalcSubROIProperties --- std::vector sys::WireModUtility::CalcSubROIProperties(sys::WireModUtility::ROIProperties_t const& roi_properties, std::vector const& hitPtrVec) { @@ -621,11 +733,12 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: void sys::WireModUtility::ModifyROI(std::vector & roi_data, sys::WireModUtility::ROIProperties_t const& roi_prop, std::vector const& subROIPropVec, - std::map const& subROIScaleMap) + std::map const& subROIScaleMap, + double sigmaWindow) { // do you want a bunch of messages? bool verbose = false; - + // initialize some values double q_orig = 0.0; double q_mod = 0.0; @@ -639,50 +752,81 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, q_mod = 0.0; scale_ratio = 1.0; + double t = i_t + roi_prop.begin; + // loop over the subs for (auto const& subroi_prop : subROIPropVec) { // get your scale vals auto scale_vals = subROIScaleMap.find(subroi_prop.key)->second; - q_orig += gausFunc(i_t + roi_prop.begin, subroi_prop.center, subroi_prop.sigma, subroi_prop.total_q); - q_mod += gausFunc(i_t + roi_prop.begin, subroi_prop.center, scale_vals.r_sigma * subroi_prop.sigma, scale_vals.r_Q * subroi_prop.total_q); + double q_this_orig = gausFunc(t, subroi_prop.center, subroi_prop.sigma, subroi_prop.total_q); + q_orig += q_this_orig; + + // sigmaWindow gating (limit scaling to within sigmaWindow of the subROI center) + // is disabled below; scale is currently applied everywhere regardless of distance. + double dist_sigma = (subroi_prop.sigma > 0) + ? std::abs(t - subroi_prop.center) / subroi_prop.sigma + : std::numeric_limits::max(); + + //if (dist_sigma <= sigmaWindow) + q_mod += gausFunc(t, subroi_prop.center, scale_vals.r_sigma * subroi_prop.sigma, scale_vals.r_Q * subroi_prop.total_q); + //else + // q_mod += q_this_orig; if (verbose) - std::cout << " Incrementing q_orig by gausFunc(" << i_t + roi_prop.begin << ", " << subroi_prop.center << ", " << subroi_prop.sigma << ", " << subroi_prop.total_q << ")" << '\n' - << " Incrementing q_mod by gausFunc(" << i_t + roi_prop.begin << ", " << subroi_prop.center << ", " << scale_vals.r_sigma * subroi_prop.sigma << ", " << scale_vals.r_Q * subroi_prop.total_q << ")" << std::endl; + std::cout << " tick " << t << " subROI center=" << subroi_prop.center + << " dist=" << dist_sigma << "sigma" + << (dist_sigma <= sigmaWindow ? " [scaled]" : " [pass-through]") << '\n' + << " q_orig+=" << q_this_orig << '\n'; } - // do some sanity checks - if (isnan(q_orig)) + if (additiveModification) { + // additive approach: after(t) = before(t) + [q_mod(t) - q_orig(t)] + // residuals/noise in the waveform are preserved without amplification + double delta = q_mod - q_orig; + if (isnan(delta) || isinf(delta)) delta = 0.0; + roi_data[i_t] += static_cast(delta); if (verbose) - std::cout << "WARNING: obtained q_orig = NaN... setting scale to 1" << std::endl; - scale_ratio = 1.0; - } else if (isnan(q_mod)) { - if (verbose) - std::cout << "WARNING: obtained q_mod = NaN... setting scale to 0" << std::endl; - scale_ratio = 0.0; - } else if (q_orig < 0.01) { // check that this is a sane limit - std::cout << "WARNING: obtained q_orig < 0.01 ... setting scale to 1" << std::endl; - scale_ratio = 1.0; - } else { - scale_ratio = q_mod / q_orig; + std::cout << "\t tick " << i_t << ":" + << " data=" << roi_data[i_t] + << ", q_orig=" << q_orig + << ", q_mod=" << q_mod + << ", delta=" << delta << " [additive]" << std::endl; } - if(isnan(scale_ratio) || isinf(scale_ratio)) + else { + // multiplicative approach (original behavior): after(t) = before(t) * [q_mod(t)/q_orig(t)] + if (isnan(q_orig)) + { + if (verbose) + std::cout << "WARNING: obtained q_orig = NaN... setting scale to 1" << std::endl; + scale_ratio = 1.0; + } else if (isnan(q_mod)) { + if (verbose) + std::cout << "WARNING: obtained q_mod = NaN... setting scale to 0" << std::endl; + scale_ratio = 0.0; + } else if (q_orig < 0.01) { + std::cout << "WARNING: obtained q_orig < 0.01 ... setting scale to 1" << std::endl; + scale_ratio = 1.0; + } else { + scale_ratio = q_mod / q_orig; + } + if(isnan(scale_ratio) || isinf(scale_ratio)) + { + if (verbose) + std::cout << "WARNING: obtained scale_ratio = " << q_mod << " / " << q_orig << " = NAN/Inf... setting to 1" << std::endl; + scale_ratio = 1.0; + } + roi_data[i_t] = scale_ratio * roi_data[i_t]; if (verbose) - std::cout << "WARNING: obtained scale_ratio = " << q_mod << " / " << q_orig << " = NAN/Inf... setting to 1" << std::endl; - scale_ratio = 1.0; + std::cout << "\t tick " << i_t << ":" + << " data=" << roi_data[i_t] + << ", q_orig=" << q_orig + << ", q_mod=" << q_mod + << ", ratio=" << scale_ratio << " [multiplicative]" << std::endl; } - - roi_data[i_t] = scale_ratio * roi_data[i_t]; - if (verbose) - std::cout << "\t tick " << i_t << ":" - << " data=" << roi_data[i_t] - << ", q_orig=" << q_orig - << ", q_mod=" << q_mod - << ", ratio=" << scale_ratio << std::endl; } // we're done now diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index e9e147269..3027980b3 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -21,6 +21,7 @@ #include "lardataobj/Simulation/SimEnergyDeposit.h" #include "lardataobj/Simulation/SimChannel.h" #include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" +#include "sbnobj/ICARUS/TPC/ChannelROI.h" #include "nusimdata/SimulationBase/MCParticle.h" #include "nusimdata/SimulationBase/MCTruth.h" @@ -44,6 +45,7 @@ namespace sys { bool applyYZAngleScale; // do we scale with YZ angle? bool applydEdXScale; // do we scale with dEdx? bool applyXXWScale; // do we scale with X vs ThXW? + bool additiveModification; // additive (true) vs multiplicative (false) ROI modification double readoutWindowTicks; // how many ticks are in the readout window? double tickOffset; // do we want an offset in the ticks? @@ -87,6 +89,7 @@ namespace sys { applyYZAngleScale(arg_ApplyYZAngleScale), applydEdXScale(arg_ApplydEdXScale), applyXXWScale(arg_ApplyXXWScale), + additiveModification(false), readoutWindowTicks(detProp.ReadOutWindowSize()), // the default A2795 (ICARUS TPC readout board) readout window is 4096 samples tickOffset(arg_TickOffset) // tick offset is for MC truth, default to zero and set only as necessary { @@ -226,6 +229,7 @@ namespace sys { // these are set in the .cc file ROIProperties_t CalcROIProperties(recob::Wire const&, size_t const&); + ROIProperties_t CalcROIProperties(recob::ChannelROI const&, size_t const&); std::vector> GetTargetROIs(sim::SimEnergyDeposit const&, double offset); std::vector> GetHitTargetROIs(recob::Hit const&); @@ -234,6 +238,9 @@ namespace sys { void FillROIMatchedHitMap(std::vector const&, std::vector const&); void FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, double offset); + void FillROIMatchedEdepMap(std::vector const&, std::vector const&, double offset); + void FillROIMatchedHitMap(std::vector const&, std::vector const&); + std::vector CalcSubROIProperties(ROIProperties_t const&, std::vector const&); std::map> MatchEdepsToSubROIs(std::vector const&, std::vector const&, double offset); @@ -247,7 +254,8 @@ namespace sys { void ModifyROI(std::vector &, ROIProperties_t const &, std::vector const&, - std::map const&); + std::map const&, + double sigmaWindow = std::numeric_limits::infinity()); }; // end class } // end namespace From 39dcea337528579b2185907234e0529a475388d0 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Tue, 30 Jun 2026 17:03:06 -0500 Subject: [PATCH 2/5] Build fixes. --- sbncode/WireMod/Utility/WireModUtility.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 3e58e329c..495354b23 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -810,7 +810,7 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, std::cout << " tick " << t << " subROI center=" << subroi_prop.center << " dist=" << sigma_distance << "sigma" << (sigma_distance <= sigmaWindow ? " [scaled]" : " [pass-through]") << '\n' - << " q_orig+=" << q_this_orig << '\n'; + << " q_orig=" << q_orig << '\n'; } double delta = q_mod - q_orig; From 643a2b1b3beaed16116bc3fad7224eb6e93c966c Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 1 Jul 2026 08:57:04 -0500 Subject: [PATCH 3/5] Add in utilities to lookup truth information from SimChannel's --- sbncode/WireMod/Utility/CMakeLists.txt | 5 +- sbncode/WireMod/Utility/WireModUtility.cc | 281 +++++++++++++++++++++- sbncode/WireMod/Utility/WireModUtility.hh | 47 +++- 3 files changed, 330 insertions(+), 3 deletions(-) diff --git a/sbncode/WireMod/Utility/CMakeLists.txt b/sbncode/WireMod/Utility/CMakeLists.txt index 3c441620f..42a65ae0d 100644 --- a/sbncode/WireMod/Utility/CMakeLists.txt +++ b/sbncode/WireMod/Utility/CMakeLists.txt @@ -17,8 +17,11 @@ cet_make_library( lardataobj::AnalysisBase lardataobj::RawData lardataobj::RecoBase + lardataobj::Simulation + nusimdata::SimulationBase + larevt::SpaceCharge sbnobj::ICARUS_TPC - lardataalg::DetectorInfo + lardataalg::DetectorInfo lardata::RecoObjects lardata::Utilities larreco::RecoAlg diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 495354b23..03ab6ece1 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -1,6 +1,10 @@ #include "sbncode/WireMod/Utility/WireModUtility.hh" #include "lardataalg/DetectorInfo/DetectorPropertiesData.h" +#include +#include +#include "TVector3.h" + //--- CalcROIProperties --- sys::WireModUtility::ROIProperties_t sys::WireModUtility::CalcROIProperties(recob::Wire const& wire, size_t const& roi_idx) { @@ -149,8 +153,105 @@ std::vector> sys::WireModUtility::GetHitTa } //--- FillROIMatchedIDEMap --- -void sys::WireModUtility::FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, double offset) +// SimChannel/IDE analogue of FillROIMatchedEdepMap. Unlike an edep, a sim::IDE already +// carries its readout channel (SimChannel::Channel()) and time (the TDCIDEMap TDC), so no +// geometric X->tick projection is needed -- the TDC is converted to a readout tick and the +// covering ROI is found directly. Builds the flat fIDEVec and the ROI->index map. +void sys::WireModUtility::FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, detinfo::DetectorClocksData const& clockData, double offset) +{ + ROIMatchedIDEMap.clear(); + fIDEVec.clear(); + + // map channel -> index into wireVec + std::unordered_map wireChannelMap; + for (size_t i_w = 0; i_w < wireVec.size(); ++i_w) + wireChannelMap[wireVec[i_w].Channel()] = i_w; + + for (auto const& simch : simchVec) + { + raw::ChannelID_t channel = simch.Channel(); + auto wcIt = wireChannelMap.find(channel); + if (wcIt == wireChannelMap.end()) + continue; + auto const& target_wire = wireVec.at(wcIt->second); + + // wire closest to the charge (used later for the SCE-undo) + geo::WireID wireID; + auto const wires = wireReadout->ChannelToWire(channel); + if (!wires.empty()) wireID = wires[0]; + + for (auto const& tdcide : simch.TDCIDEMap()) + { + // convert the SimChannel TDC to a readout tick (same frame as the ROI samples) + double tick = clockData.TPCTDC2Tick(tdcide.first) + offset + tickOffset; + if (tick < 0) continue; + unsigned int sample = (unsigned int) std::round(tick); + + if (not target_wire.SignalROI().is_valid() || + target_wire.SignalROI().empty() || + target_wire.SignalROI().n_ranges() == 0 || + target_wire.SignalROI().size() <= sample || + target_wire.SignalROI().is_void(sample) ) + continue; + + auto range_number = target_wire.SignalROI().find_range_iterator(sample) - target_wire.SignalROI().begin_range(); + ROI_Key_t roi_key = std::make_pair(target_wire.Channel(), range_number); + + for (sim::IDE const& ide : tdcide.second) + { + fIDEVec.push_back( MatchedIDE_t{channel, wireID, tick, &ide} ); + ROIMatchedIDEMap[roi_key].push_back(fIDEVec.size() - 1); + } + } + } +} + +//--- FillROIMatchedIDEMap (ChannelROI overload) --- +void sys::WireModUtility::FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& chanROIVec, detinfo::DetectorClocksData const& clockData, double offset) { + ROIMatchedIDEMap.clear(); + fIDEVec.clear(); + + std::unordered_map chanROIChannelMap; + for (size_t i_c = 0; i_c < chanROIVec.size(); ++i_c) + chanROIChannelMap[chanROIVec[i_c].Channel()] = i_c; + + for (auto const& simch : simchVec) + { + raw::ChannelID_t channel = simch.Channel(); + auto ccIt = chanROIChannelMap.find(channel); + if (ccIt == chanROIChannelMap.end()) + continue; + auto const& target_chanROI = chanROIVec.at(ccIt->second); + auto const& roiShort = target_chanROI.SignalROI(); + + geo::WireID wireID; + auto const wires = wireReadout->ChannelToWire(channel); + if (!wires.empty()) wireID = wires[0]; + + for (auto const& tdcide : simch.TDCIDEMap()) + { + double tick = clockData.TPCTDC2Tick(tdcide.first) + offset + tickOffset; + if (tick < 0) continue; + unsigned int sample = (unsigned int) std::round(tick); + + if (not roiShort.is_valid() || + roiShort.empty() || + roiShort.n_ranges() == 0 || + roiShort.size() <= sample || + roiShort.is_void(sample) ) + continue; + + auto range_number = roiShort.find_range_iterator(sample) - roiShort.begin_range(); + ROI_Key_t roi_key = std::make_pair(target_chanROI.Channel(), range_number); + + for (sim::IDE const& ide : tdcide.second) + { + fIDEVec.push_back( MatchedIDE_t{channel, wireID, tick, &ide} ); + ROIMatchedIDEMap[roi_key].push_back(fIDEVec.size() - 1); + } + } + } } //--- FillROIMatchedEdepMap --- @@ -600,6 +701,184 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd return edep_col_properties; } +//--- WireToTrajectoryPosition (SCE: at-the-wire -> true trajectory frame) --- +geo::Point_t sys::WireModUtility::WireToTrajectoryPosition(const geo::Point_t& loc, const geo::TPCID& tpc) const +{ + geo::Point_t ret = loc; + if (fSCE && fSCE->EnableSimSpatialSCE()) + { + geo::Vector_t offset = fSCE->GetCalPosOffsets(ret, tpc.TPC); + ret.SetX(ret.X() + offset.X()); + ret.SetY(ret.Y() + offset.Y()); + ret.SetZ(ret.Z() + offset.Z()); + } + return ret; +} + +//--- TrajectoryToWirePosition (SCE: true -> at-the-wire frame) --- +geo::Point_t sys::WireModUtility::TrajectoryToWirePosition(const geo::Point_t& loc, const geo::Vector_t& driftdir) const +{ + geo::Point_t ret = loc; + int corr = driftdir.X(); // returned X offset is the drift; sign it by the drift direction + if (fSCE && fSCE->EnableSimSpatialSCE()) + { + geo::Vector_t offset = fSCE->GetPosOffsets(ret); + ret.SetX(ret.X() + corr * offset.X()); + ret.SetY(ret.Y() + offset.Y()); + ret.SetZ(ret.Z() + offset.Z()); + } + return ret; +} + +//--- MatchIDEsToSubROIs --- +// SimChannel/IDE analogue of MatchEdepsToSubROIs. Each IDE already carries its readout tick +// (no per-plane X->tick projection needed), so it is assigned to the closest sub-ROI whose +// distance |tick - center|/sigma is below the same threshold used for edeps (5 sigma). +std::map> +sys::WireModUtility::MatchIDEsToSubROIs(std::vector const& subROIPropVec, + std::vector const& idePtrVec) +{ + std::map> SubROIMatchedIDEMap; // subROI idx -> ide idxs + + for (unsigned int i_e = 0; i_e < idePtrVec.size(); ++i_e) + { + double ide_tick = idePtrVec[i_e]->tick; + unsigned int closest = std::numeric_limits::max(); + float min_dist = std::numeric_limits::max(); + for (unsigned int i_h = 0; i_h < subROIPropVec.size(); ++i_h) + { + auto const& s = subROIPropVec[i_h]; + if (s.sigma <= 0) continue; + float dist = std::abs((float)(ide_tick - s.center)) / s.sigma; + if (dist < min_dist) { min_dist = dist; closest = i_h; } + } + if (closest != std::numeric_limits::max() && min_dist < 5) + SubROIMatchedIDEMap[closest].push_back(i_e); + } + + std::map> ReturnMap; + for (auto const& kv : SubROIMatchedIDEMap) + { + auto key = subROIPropVec[kv.first].key; + for (auto const& i_e : kv.second) + ReturnMap[key].push_back(idePtrVec[i_e]); + } + return ReturnMap; +} + +//--- CalcPropertiesFromIDEs --- +// SimChannel/IDE analogue of CalcPropertiesFromEdeps. Picks the dominant track (max summed +// energy, abs(trackID) as in TrackCaloSkimmer), builds the numElectrons-weighted centroid in +// the SCE-undone (trajectory) frame, and reads the direction (dxdr/dydr/dzdr) from the matched +// MCParticle's nearest trajectory point. dedr/dqdr use the trajectory-derived pitch. The +// scales_avg pre-pass of the edep version is dropped (those fields are never read downstream). +sys::WireModUtility::TruthProperties_t +sys::WireModUtility::CalcPropertiesFromIDEs(std::vector const& idePtrVec, + std::map const& particleMap) +{ + TruthProperties_t props; + props.x = props.y = props.z = 0.f; + props.x_rms = props.x_rms_noWeight = 0.f; + props.x_min = std::numeric_limits::max(); + props.x_max = std::numeric_limits::lowest(); + props.tick = props.tick_rms = props.tick_rms_noWeight = 0.f; + props.tick_min = std::numeric_limits::max(); + props.tick_max = std::numeric_limits::lowest(); + props.dxdr = props.dydr = props.dzdr = 0.; + props.dedr = props.dqdr = 0.; + props.total_energy = 0.f; + for (size_t i_p = 0; i_p < 3; ++i_p) { props.scales_avg[i_p].r_Q = 1.; props.scales_avg[i_p].r_sigma = 1.; } + if (idePtrVec.empty()) return props; + + // dominant track = max summed energy (abs(trackID), matching TrackCaloSkimmer) + std::map energy_per_trkid; + for (auto const& m : idePtrVec) + energy_per_trkid[std::abs(m->ide->trackID)] += m->ide->energy; + int trkid_max = 0; double energy_max = -1.; + for (auto const& e : energy_per_trkid) + if (e.second > energy_max) { energy_max = e.second; trkid_max = e.first; } + + // numElectrons-weighted centroid in the SCE-undone frame + energy/electron/tick sums + struct Acc { double xtraj, tick, ne; }; + std::vector accs; + double sumNe = 0., sumE = 0., cx = 0., cy = 0., cz = 0., tsum = 0.; + geo::WireID domWire; + for (auto const& m : idePtrVec) + { + if (std::abs(m->ide->trackID) != trkid_max) continue; + double ne = m->ide->numElectrons; + double e = m->ide->energy; + geo::Point_t p_raw(m->ide->x, m->ide->y, m->ide->z); + geo::Point_t p_traj = WireToTrajectoryPosition(p_raw, m->wire); + cx += ne * p_traj.X(); cy += ne * p_traj.Y(); cz += ne * p_traj.Z(); + tsum += ne * m->tick; + sumNe += ne; sumE += e; + accs.push_back({p_traj.X(), m->tick, ne}); + if (p_traj.X() < props.x_min) props.x_min = p_traj.X(); + if (p_traj.X() > props.x_max) props.x_max = p_traj.X(); + if (m->tick < props.tick_min) props.tick_min = m->tick; + if (m->tick > props.tick_max) props.tick_max = m->tick; + domWire = m->wire; + } + if (sumNe <= 0) return props; + cx /= sumNe; cy /= sumNe; cz /= sumNe; + props.x = cx; props.y = cy; props.z = cz; + props.total_energy = sumE; + props.tick = tsum / sumNe; + + double xr = 0., xrnw = 0., tr = 0.; + for (auto const& a : accs) + { + xr += a.ne * (a.xtraj - cx) * (a.xtraj - cx); + xrnw += (a.xtraj - cx) * (a.xtraj - cx); + tr += a.ne * (a.tick - props.tick) * (a.tick - props.tick); + } + props.x_rms = std::sqrt(xr / sumNe); + props.x_rms_noWeight = std::sqrt(xrnw); + props.tick_rms = std::sqrt(tr / sumNe); + props.tick_rms_noWeight = props.x_rms_noWeight; + + // direction / pitch from the matched MCParticle trajectory (nearest point to centroid) + auto pit = particleMap.find(trkid_max); + if (pit != particleMap.end() && pit->second != nullptr) + { + const simb::MCParticle* part = pit->second; + TVector3 hp(cx, cy, cz); + double closest = -1.; TVector3 dir; + for (unsigned int i = 0; i < part->NumberTrajectoryPoints(); ++i) + { + double d = (part->Position(i).Vect() - hp).Mag(); + if (closest < 0. || d < closest) + { + closest = d; + dir = part->Momentum(i).Vect().Unit(); + } + } + if (closest >= 0. && dir.Mag() > 1e-4) + { + props.dxdr = dir.X(); + props.dydr = dir.Y(); + props.dzdr = dir.Z(); + + if (domWire.isValid) + { + geo::PlaneID plane(domWire.Cryostat, domWire.TPC, domWire.Plane); + geo::PlaneGeo const& pg = wireReadout->Plane(plane); + double angletovert = wireReadout->WireAngleToVertical(pg.View(), plane) - 0.5 * ::util::pi<>(); + double cosgamma = std::abs(std::cos(angletovert) * dir.Z() + std::sin(angletovert) * dir.Y()); + double pitch = (cosgamma > 1e-6) ? pg.WirePitch() / cosgamma : pg.WirePitch(); + if (pitch > 0) + { + props.dedr = sumE / pitch; // dE/dx [MeV/cm] + props.dqdr = sumNe / pitch; // dQ/dx [electrons/cm] + } + } + } + } + + return props; +} + //--- GetScaleValues --- sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetScaleValues(sys::WireModUtility::TruthProperties_t const& truth_props, sys::WireModUtility::ROIProperties_t const& roi_vals) { diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index ba1d15fa0..6ca76c9e3 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -14,6 +14,8 @@ #include "larcorealg/Geometry/GeometryCore.h" #include "larcorealg/Geometry/WireReadoutGeom.h" #include "lardataalg/DetectorInfo/DetectorPropertiesData.h" +#include "lardataalg/DetectorInfo/DetectorClocksData.h" +#include "larevt/SpaceCharge/SpaceCharge.h" #include "lardataobj/RecoBase/Track.h" #include "lardataobj/RecoBase/TrackHitMeta.h" #include "lardataobj/RecoBase/Hit.h" @@ -166,10 +168,32 @@ namespace sys { ScaleValues_t scales_avg[3]; } TruthProperties_t; + // A single SimChannel ionization deposit (sim::IDE) tagged with its readout + // coordinates. This is the SimChannel analogue of a sim::SimEnergyDeposit used by + // the truth-matching path; unlike an edep, a sim::IDE carries only a point position + // (x,y,z) -- direction/step come from the matched simb::MCParticle trajectory. + typedef struct MatchedIDE + { + raw::ChannelID_t channel; // readout channel the charge landed on + geo::WireID wire; // wire closest to the deposit + double tick; // readout tick (TDC converted to tick) + const sim::IDE* ide; // trackID, numElectrons, energy, x, y, z + } MatchedIDE_t; + // internal containers std::map< ROI_Key_t, std::vector > ROIMatchedEdepMap; std::map< ROI_Key_t, std::vector > ROIMatchedHitMap; + // SimChannel/IDE truth-matching containers (parallel to the Edep versions). + // fIDEVec is the flat owner; ROIMatchedIDEMap stores indices into it, keyed by ROI. + std::vector fIDEVec; + std::map< ROI_Key_t, std::vector > ROIMatchedIDEMap; + + // Space-charge provider, set by the module (lar::providerFrom()). + // Used to map IDE (at-the-wire) positions to/from the true MCParticle trajectory + // frame. When null or SCE disabled, the mappings are the identity. + const spacecharge::SpaceCharge* fSCE = nullptr; + // some useful functions // geometries // TODO is this the most efficient for new v10 iterators? @@ -254,7 +278,12 @@ namespace sys { void FillROIMatchedEdepMap(std::vector const&, std::vector const&, double offset); void FillROIMatchedHitMap(std::vector const&, std::vector const&); - void FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, double offset); + + // SimChannel/IDE versions of FillROIMatchedEdepMap. The DetectorClocksData is used to + // convert each IDE's TDC to a readout tick (no X->tick projection / tickOffset needed, + // since the SimChannel already carries the readout coordinate). + void FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, detinfo::DetectorClocksData const& clockData, double offset = 0); + void FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& chanROIVec, detinfo::DetectorClocksData const& clockData, double offset = 0); void FillROIMatchedEdepMap(std::vector const&, std::vector const&, double offset); void FillROIMatchedHitMap(std::vector const&, std::vector const&); @@ -263,8 +292,24 @@ namespace sys { std::map> MatchEdepsToSubROIs(std::vector const&, std::vector const&, double offset); + // SimChannel/IDE analogue of MatchEdepsToSubROIs. Keyed on the IDE's native readout + // tick (no per-plane X->tick projection), otherwise the same center+/-sigma / closest + // sub-ROI assignment as the edep version. + std::map> MatchIDEsToSubROIs(std::vector const&, std::vector const&); + TruthProperties_t CalcPropertiesFromEdeps(std::vector const&, double offset); + // SimChannel/IDE analogue of CalcPropertiesFromEdeps. Fills TruthProperties_t from the + // dominant track's sim::IDEs (charge-weighted position/energy) plus the matched + // simb::MCParticle trajectory (direction / pitch / dedr / dqdr). particleMap maps + // abs(trackID) -> MCParticle so the dominant track's trajectory can be looked up. + TruthProperties_t CalcPropertiesFromIDEs(std::vector const&, std::map const& particleMap); + + // Space-charge coordinate maps (mirroring TrackCaloSkimmer). Identity if fSCE is null + // or SCE disabled. + geo::Point_t WireToTrajectoryPosition(const geo::Point_t& loc, const geo::TPCID& tpc) const; // at-the-wire -> true trajectory frame + geo::Point_t TrajectoryToWirePosition(const geo::Point_t& loc, const geo::Vector_t& driftdir) const; // true -> at-the-wire frame + ScaleValues_t GetScaleValues(TruthProperties_t const&, ROIProperties_t const&); ScaleValues_t GetChannelScaleValues(TruthProperties_t const&, raw::ChannelID_t const&); ScaleValues_t GetViewScaleValues(TruthProperties_t const&, geo::View_t const&); From b110852e869ff7f9db0ae5dc199446636b164da5 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Thu, 2 Jul 2026 14:08:44 -0500 Subject: [PATCH 4/5] Fix thXW calculation. --- sbncode/WireMod/Utility/WireModUtility.hh | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index 6ca76c9e3..3755a1b0a 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -260,9 +260,14 @@ namespace sys { double ThetaXW(double dxdr, double dydr, double dzdr, double planeAngle) { - // Want angle to vert, so subtract pi - double sin = std::sin(planeAngle - ::util::pi<>()); - double cos = std::cos(planeAngle - ::util::pi<>()); + // planeAngle is the wire angle from +z (PlaneGeo::ThetaZ, numerically equal to + // WireReadoutGeom::WireAngleToVertical). The pitch "gamma" angle is measured from the + // wire-normal, i.e. planeAngle - pi/2 -- the same convention as the canonical LArSoft + // pitch (Calorimetry / TrackCaloSkimmer use WireAngleToVertical - pi/2). Subtracting a + // full pi projected the direction onto the wire itself (a cos<->sin swap), yielding the + // wrong angle; subtract pi/2 so cosG is the projection onto the wire-normal (pitch). + double sin = std::sin(planeAngle - 0.5 * ::util::pi<>()); + double cos = std::cos(planeAngle - 0.5 * ::util::pi<>()); double cosG = std::abs(dydr * sin + dzdr * cos); double theta = std::atan(dxdr / cosG); From 2b53154f7750221ed5be67df00ec1ff4857265f6 Mon Sep 17 00:00:00 2001 From: "poppi@baltig.infn.it" Date: Wed, 8 Jul 2026 23:53:59 +0200 Subject: [PATCH 5/5] Some fixes to wiremod utility. First we removed the modification cutoff when using additive modification. Second I modified the interpolation that was stopping at bin center for bins at the borders, in this case we clamp the interpolation to the outer bin edge. Third I introduced a flag which allows WireMod to either use interpolation or pick up the scaling factor from the closest bin center. --- sbncode/WireMod/Utility/WireModUtility.cc | 76 ++++++++++++++++++----- sbncode/WireMod/Utility/WireModUtility.hh | 2 + 2 files changed, 61 insertions(+), 17 deletions(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 03ab6ece1..001ef7ba0 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -1,6 +1,7 @@ #include "sbncode/WireMod/Utility/WireModUtility.hh" #include "lardataalg/DetectorInfo/DetectorPropertiesData.h" +#include #include #include #include "TVector3.h" @@ -879,6 +880,49 @@ sys::WireModUtility::CalcPropertiesFromIDEs(std::vectorInterpolate(std::clamp(x, g->GetXmin(), g->GetXmax()), + std::clamp(y, g->GetYmin(), g->GetYmax())); + + int n = g->GetN(); + double* xs = g->GetX(); + double* ys = g->GetY(); + double* zs = g->GetZ(); + + std::vector ux(xs, xs + n), uy(ys, ys + n); + std::sort(ux.begin(), ux.end()); + ux.erase(std::unique(ux.begin(), ux.end()), ux.end()); + std::sort(uy.begin(), uy.end()); + uy.erase(std::unique(uy.begin(), uy.end()), uy.end()); + + double qx = std::clamp(x, ux.front(), ux.back()); + double qy = std::clamp(y, uy.front(), uy.back()); + + auto findBin = [](std::vector const& c, double q) -> int { + int m = static_cast(c.size()); + if (m == 1) return 0; + double left = c[0] - (c[1] - c[0]) * 0.5; + for (int i = 0; i < m; ++i) { + double right = 2.0 * c[i] - left; + if (q < right || i == m - 1) return i; + left = right; + } + return m - 1; + }; + + double cx = ux[findBin(ux, qx)]; + double cy = uy[findBin(uy, qy)]; + + for (int i = 0; i < n; ++i) + if (std::abs(xs[i] - cx) < 1e-9 && std::abs(ys[i] - cy) < 1e-9) + return zs[i]; + + return 1.0; // bin not populated in the map (sparse grid) — neutral scale +} + //--- GetScaleValues --- sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetScaleValues(sys::WireModUtility::TruthProperties_t const& truth_props, sys::WireModUtility::ROIProperties_t const& roi_vals) { @@ -955,10 +999,10 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: graph2Ds_Sigma_YZ [plane] == nullptr ) throw cet::exception("WireModUtility") << "Tried to apply YZ scale factor, but could not find graphs. Check that you have set those in the utility."; - temp_scale = graph2Ds_Charge_YZ[plane]->Interpolate(truth_props.z, truth_props.y); + temp_scale = g2dLookup(graph2Ds_Charge_YZ[plane], truth_props.z, truth_props.y, useGraph2DInterpolation); if(temp_scale>0.001) scales.r_Q *= temp_scale; - temp_scale = graph2Ds_Sigma_YZ [plane]->Interpolate(truth_props.z, truth_props.y); + temp_scale = g2dLookup(graph2Ds_Sigma_YZ [plane], truth_props.z, truth_props.y, useGraph2DInterpolation); if(temp_scale>0.001) scales.r_sigma *= temp_scale; } @@ -969,10 +1013,10 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: throw cet::exception("WireModUtility") << "Tried to apply XXW scale factor, but could not find graphs. Check that you have set those in the utility."; double thXW = ThetaXW(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()); - temp_scale = graph2Ds_Charge_XXW[plane]->Interpolate(truth_props.x, thXW); + temp_scale = g2dLookup(graph2Ds_Charge_XXW[plane], truth_props.x, thXW, useGraph2DInterpolation); if(temp_scale>0.001) scales.r_Q *= temp_scale; - temp_scale = graph2Ds_Sigma_XXW [plane]->Interpolate(truth_props.x, thXW); + temp_scale = g2dLookup(graph2Ds_Sigma_XXW [plane], truth_props.x, thXW, useGraph2DInterpolation); if(temp_scale>0.001) scales.r_sigma *= temp_scale; } @@ -1011,10 +1055,10 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: graph2Ds_Sigma_XXZAngle [plane] == nullptr ) throw cet::exception("WireModUtility") << "Tried to apply XXZAngle scale factor, but could not find graphs. Check that you have set those in the utility."; - temp_scale = graph2Ds_Charge_XXZAngle[plane]->Interpolate(truth_props.x, ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + temp_scale = g2dLookup(graph2Ds_Charge_XXZAngle[plane], truth_props.x, ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), useGraph2DInterpolation); if(temp_scale>0.001) scales.r_Q *= temp_scale; - temp_scale = graph2Ds_Sigma_XXZAngle [plane]->Interpolate(truth_props.x, ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + temp_scale = g2dLookup(graph2Ds_Sigma_XXZAngle [plane], truth_props.x, ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), useGraph2DInterpolation); if(temp_scale>0.001) scales.r_sigma *= temp_scale; } @@ -1024,10 +1068,10 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: graph2Ds_Sigma_XdQdX [plane] == nullptr ) throw cet::exception("WireModUtility") << "Tried to apply XdQdX scale factor, but could not find graphs. Check that you have set those in the utility."; - temp_scale = graph2Ds_Charge_XdQdX[plane]->Interpolate(truth_props.x, truth_props.dqdr); + temp_scale = g2dLookup(graph2Ds_Charge_XdQdX[plane], truth_props.x, truth_props.dqdr, useGraph2DInterpolation); if(temp_scale>0.001) scales.r_Q *= temp_scale; - temp_scale = graph2Ds_Sigma_XdQdX [plane]->Interpolate(truth_props.x, truth_props.dqdr); + temp_scale = g2dLookup(graph2Ds_Sigma_XdQdX [plane], truth_props.x, truth_props.dqdr, useGraph2DInterpolation); if(temp_scale>0.001) scales.r_sigma *= temp_scale; } @@ -1037,10 +1081,10 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: graph2Ds_Sigma_XZAngledQdX [plane] == nullptr ) throw cet::exception("WireModUtility") << "Tried to apply XZAngledQdX scale factor, but could not find graphs. Check that you have set those in the utility."; - temp_scale = graph2Ds_Charge_XZAngledQdX[plane]->Interpolate(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), truth_props.dqdr); + temp_scale = g2dLookup(graph2Ds_Charge_XZAngledQdX[plane], ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), truth_props.dqdr, useGraph2DInterpolation); if(temp_scale>0.001) scales.r_Q *= temp_scale; - temp_scale = graph2Ds_Sigma_XZAngledQdX [plane]->Interpolate(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), truth_props.dqdr); + temp_scale = g2dLookup(graph2Ds_Sigma_XZAngledQdX [plane], ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), truth_props.dqdr, useGraph2DInterpolation); if(temp_scale>0.001) scales.r_sigma *= temp_scale; } @@ -1098,16 +1142,14 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, if (isnan(q_mod)) { if (verbose) std::cout << "WARNING: obtained q_mod = NaN... setting scale to 0" << std::endl; + // for additive correction, no cutoffs are applied + } else if (additiveModification) { + roi_data[i_t] += static_cast(delta); } else if (q_orig < 0.01) { // check that this is a sane limit if (verbose) std::cout << "WARNING: obtained q_orig < 0.01 ... setting scale to 1" << std::endl; } else if (sigma_distance > 9.) { - if (verbose) std::cout << "WARNING: sigma_distance > 3.. setting scale to 1" << std::endl; - // additive approach: after(t) = before(t) + [q_mod(t) - q_orig(t)] - // residuals/noise in the waveform are preserved without amplification - } else if (additiveModification) { - roi_data[i_t] += static_cast(delta); - } - else { + if (verbose) std::cout << "WARNING: sigma_distance > 9 ... setting scale to 1" << std::endl; + } else { scale_ratio = q_mod / q_orig; roi_data[i_t] = scale_ratio * roi_data[i_t]; } diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index 3755a1b0a..885a248f0 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -51,6 +51,7 @@ namespace sys { bool applyXZAngledQdXScale; // do we scale with XZ angle vs dQ/dX? bool applyXXWScale; // do we scale with X vs ThXW? bool additiveModification; // additive (true) vs multiplicative (false) ROI modification + bool useGraph2DInterpolation; // true = Delaunay interp (default), false = nearest bin center double readoutWindowTicks; // how many ticks are in the readout window? double tickOffset; // do we want an offset in the ticks? @@ -108,6 +109,7 @@ namespace sys { applyXZAngledQdXScale(arg_ApplyXZAngledQdXScale), applyXXWScale(arg_ApplyXXWScale), additiveModification(false), + useGraph2DInterpolation(true), readoutWindowTicks(detProp.ReadOutWindowSize()), // the default A2795 (ICARUS TPC readout board) readout window is 4096 samples tickOffset(arg_TickOffset) // tick offset is for MC truth, default to zero and set only as necessary {