From 41e4fa4d0a203b456479b3f0d605c99bf55306b6 Mon Sep 17 00:00:00 2001 From: Harry Hausner Date: Wed, 5 Jun 2024 14:10:18 -0500 Subject: [PATCH 01/21] add wiremod to sbncode --- sbncode/CMakeLists.txt | 1 + sbncode/WireMod/CMakeLists.txt | 1 + sbncode/WireMod/Utility/CMakeLists.txt | 45 ++ sbncode/WireMod/Utility/WireModUtility.cc | 632 ++++++++++++++++++++++ sbncode/WireMod/Utility/WireModUtility.hh | 218 ++++++++ 5 files changed, 897 insertions(+) create mode 100644 sbncode/WireMod/CMakeLists.txt create mode 100644 sbncode/WireMod/Utility/CMakeLists.txt create mode 100644 sbncode/WireMod/Utility/WireModUtility.cc create mode 100644 sbncode/WireMod/Utility/WireModUtility.hh diff --git a/sbncode/CMakeLists.txt b/sbncode/CMakeLists.txt index daed2590a..8b9a8ff9f 100644 --- a/sbncode/CMakeLists.txt +++ b/sbncode/CMakeLists.txt @@ -21,6 +21,7 @@ add_subdirectory(CosmicID) add_subdirectory(DetSim) add_subdirectory(Cluster3D) add_subdirectory(HitFinder) +add_subdirectory(WireMod) # Supera # diff --git a/sbncode/WireMod/CMakeLists.txt b/sbncode/WireMod/CMakeLists.txt new file mode 100644 index 000000000..9919ce928 --- /dev/null +++ b/sbncode/WireMod/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(Utility) diff --git a/sbncode/WireMod/Utility/CMakeLists.txt b/sbncode/WireMod/Utility/CMakeLists.txt new file mode 100644 index 000000000..8955de261 --- /dev/null +++ b/sbncode/WireMod/Utility/CMakeLists.txt @@ -0,0 +1,45 @@ +cet_make_library( + LIBRARY_NAME + sbncode_WireMod_Utility + + SOURCE + WireModUtility.cc + + LIBRARIES + ${ART_FRAMEWORK_CORE} + ${MF_MESSAGELOGGER} + ${Boost_SYSTEM_LIBRARY} + ${ROOT_BASIC_LIB_LIST} + ${ROOT_HIST} + art_root_io::TFileService_service + larcorealg::Geometry + larcore::Geometry_Geometry_service + lardataobj::AnalysisBase + lardataobj::RawData + lardataobj::RecoBase + lardataalg::DetectorInfo + lardata::RecoObjects + lardata::Utilities + larreco::RecoAlg + ${LARRECO_LIB} + ${LARDATA_LIB} + ${ART_FRAMEWORK_CORE} + ${ART_FRAMEWORK_PRINCIPAL} + ${ART_FRAMEWORK_SERVICES_REGISTRY} + ${ART_FRAMEWORK_SERVICES_OPTIONAL} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_RANDOMNUMBERGENERATOR_SERVICE} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_TFILESERVICE_SERVICE} + ${MF_MESSAGELOGGER} + ${FHICLCPP} + ${CLHEP} + ${ROOT_GEOM} + ${ROOT_XMLIO} + ${ROOT_GDML} + ${ROOT_BASIC_LIB_LIST} + ${ROOT_HIST} + BASENAME_ONLY +) + +install_headers() +install_source() + diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc new file mode 100644 index 000000000..71b0bb3ed --- /dev/null +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -0,0 +1,632 @@ +#include "sbncode/WireMod/Utility/WireModUtility.hh" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" + +//--- CalcROIProperties --- +sys::WireModUtility::ROIProperties_t sys::WireModUtility::CalcROIProperties(recob::Wire const& wire, size_t const& roi_idx) +{ + // get the ROI + recob::Wire::RegionsOfInterest_t::datarange_t const& roi = wire.SignalROI().get_ranges()[roi_idx]; + + // initialize the return value + ROIProperties_t roi_vals; + roi_vals.channel = wire.Channel(); + roi_vals.view = wire.View(); + roi_vals.begin = roi.begin_index(); + roi_vals.end = roi.end_index(); + roi_vals.center = 0; + roi_vals.total_q = 0; + roi_vals.sigma = 0; + + // loop over the roi and find the charge-weighted center and total charge + auto const& roi_data = roi.data(); + for (size_t i_t = 0; i_t < roi_data.size(); ++i_t) + { + roi_vals.center += roi_data[i_t]*(i_t+roi_vals.begin); + roi_vals.total_q += roi_data[i_t]; + } + roi_vals.center = roi_vals.center/roi_vals.total_q; + + // get the width (again charge-weighted) + // if the ROI is only one tick set the cent to the middle of the tick and width to 0.5 + for (size_t i_t = 0; i_t> sys::WireModUtility::GetTargetROIs(sim::SimEnergyDeposit const& shifted_edep, double offset) +{ + // initialize return value + // the pairs are + std::vector> target_roi_vec; + + // if the SimEnergyDeposit isn't inside the TPC then it's not going to match a wire + geo::TPCGeo const* curTPCGeomPtr = geometry->PositionToTPCptr(shifted_edep.MidPoint()); + if (curTPCGeomPtr == nullptr) + return target_roi_vec; + + // iterate over planes + for (size_t i_p = 0; i_p < curTPCGeomPtr->Nplanes(); ++i_p) + { + // check the wire exists in this plane + int wireNumber = int(0.5 + curTPCGeomPtr->Plane(i_p).WireCoordinate(shifted_edep.MidPoint())); + if ((wireNumber < 0) || (wireNumber >= (int) curTPCGeomPtr->Plane(i_p).Nwires())) + continue; + + // reconstruct the wireID from the position + geo::WireID edep_wireID = curTPCGeomPtr->Plane(i_p).NearestWireID(shifted_edep.MidPoint()); + + // if the deposition is inside the range where it would leave a signal on the wire, look for the ROI there + if (planeXInWindow(shifted_edep.X(), i_p, *curTPCGeomPtr, offset + tickOffset)) + target_roi_vec.emplace_back(geometry->PlaneWireToChannel(edep_wireID), std::round(planeXToTick(shifted_edep.X(), i_p, *curTPCGeomPtr, offset + tickOffset))); + } + + return target_roi_vec; +} + +//--- GetHitTargetROIs --- +std::vector> sys::WireModUtility::GetHitTargetROIs(recob::Hit const& hit) +{ + // initialize return value + // the pairs are + std::vector> target_roi_vec; + + int hit_wire = hit.Channel(); + int hit_tick = int(round(hit.PeakTime())); + + if (hit_tick < tickOffset || hit_tick >= readoutWindowTicks + tickOffset) + return target_roi_vec; + + + target_roi_vec.emplace_back((unsigned int) hit_wire, (unsigned int) hit_tick); + + return target_roi_vec; +} + +//--- FillROIMatchedEdepMap --- +void sys::WireModUtility::FillROIMatchedEdepMap(std::vector const& edepVec, std::vector const& wireVec, double offset) +{ + // clear the map in case it was already set + ROIMatchedEdepMap.clear(); + + // get the channel from each wire and set up a map between them + std::unordered_map wireChannelMap; + for (size_t i_w = 0; i_w < wireVec.size(); ++i_w) + wireChannelMap[wireVec[i_w].Channel()] = i_w; + + // loop over the energy deposits + for (size_t i_e = 0; i_e < edepVec.size(); ++i_e) + { + // get the ROIs + // + std::vector> target_rois = GetTargetROIs(edepVec[i_e], offset); + + // loop over ROI and match the energy deposits with wires + for (auto const& target_roi : target_rois) + { + // if we can't find the wire, skip it + if (wireChannelMap.find(target_roi.first) == wireChannelMap.end()) + continue; + + // get the wire + auto const& target_wire = wireVec.at(wireChannelMap[target_roi.first]); + + // if there are no ticks in in the wire skip it + // likewise if there's nothing in the region of interst + if (not target_wire.SignalROI().is_valid() || + target_wire.SignalROI().empty() || + target_wire.SignalROI().n_ranges() == 0 || + target_wire.SignalROI().size() <= target_roi.second || + target_wire.SignalROI().is_void(target_roi.second) ) + continue; + + // how far into the range is the ROI? + auto range_number = target_wire.SignalROI().find_range_iterator(target_roi.second) - target_wire.SignalROI().begin_range(); + + // popluate the map + ROIMatchedEdepMap[std::make_pair(target_wire.Channel(),range_number)].push_back(i_e); + } + } +} + +//--- FillROIMatchedHitMap --- +void sys::WireModUtility::FillROIMatchedHitMap(std::vector const& hitVec, std::vector const& wireVec) +{ + // clear the map in case it was already set + ROIMatchedHitMap.clear(); + + // get the channel from each wire and set up a map between them + std::unordered_map wireChannelMap; + for (size_t i_w = 0; i_w < wireVec.size(); ++i_w) + wireChannelMap[wireVec[i_w].Channel()] = i_w; + + // loop over hits + for (size_t i_h = 0; i_h < hitVec.size(); ++i_h) + { + // get the ROIs + // + std::vector> target_rois = GetHitTargetROIs(hitVec[i_h]); + + // loop over ROI and match the energy deposits with wires + for (auto const& target_roi : target_rois) + { + // if we can't find the wire, skip it + if (wireChannelMap.find(target_roi.first) == wireChannelMap.end()) + continue; + + auto const& target_wire = wireVec.at(wireChannelMap[target_roi.first]); + + // if there are no ticks in in the wire skip it + // likewise if there's nothing in the region of interst + if (not target_wire.SignalROI().is_valid() || + target_wire.SignalROI().empty() || + target_wire.SignalROI().n_ranges() == 0 || + target_wire.SignalROI().size() <= target_roi.second || + target_wire.SignalROI().is_void(target_roi.second) ) + continue; + + // which range is it? + auto range_number = target_wire.SignalROI().find_range_iterator(target_roi.second) - target_wire.SignalROI().begin_range(); + + // pupluate the map + ROIMatchedHitMap[std::make_pair(target_wire.Channel(),range_number)].push_back(i_h); + } + } +} + +//--- CalcSubROIProperties --- +std::vector sys::WireModUtility::CalcSubROIProperties(sys::WireModUtility::ROIProperties_t const& roi_properties, std::vector const& hitPtrVec) +{ + std::vector subroi_properties_vec; + sys::WireModUtility::SubROIProperties_t subroi_properties; + subroi_properties.channel = roi_properties.channel; + subroi_properties.view = roi_properties.view; + + // if this ROI doesn't contain any hits, define subROI based on ROI properities + // otherwise, define subROIs based on hits + if (hitPtrVec.size() == 0) + { + subroi_properties.key = std::make_pair(roi_properties.key, 0); + subroi_properties.total_q = roi_properties.total_q; + subroi_properties.center = roi_properties.center; + subroi_properties.sigma = roi_properties.sigma; + subroi_properties_vec.push_back(subroi_properties); + } else + { + for (unsigned int i_h=0; i_h < hitPtrVec.size(); ++i_h) + { + auto hit_ptr = hitPtrVec[i_h]; + subroi_properties.key = std::make_pair(roi_properties.key, i_h); + subroi_properties.total_q = hit_ptr->Integral(); + subroi_properties.center = hit_ptr->PeakTime(); + subroi_properties.sigma = hit_ptr->RMS(); + subroi_properties_vec.push_back(subroi_properties); + } + } + + return subroi_properties_vec; +} + +//--- MatchEdepsToSubROIs --- +std::map> sys::WireModUtility::MatchEdepsToSubROIs(std::vector const& subROIPropVec, + std::vector const& edepPtrVec, double offset) +{ + // for each TrackID, which EDeps are associated with it? keys are TrackIDs + std::map> TrackIDMatchedEDepMap; + + // total energy of EDeps matched to the ROI (not strictly necessary, but useful for understanding/development + double total_energy = 0.0; + + // loop over edeps, fill TrackIDMatchedEDepMap and calculate total energy + for (auto const& edep_ptr : edepPtrVec) + { + TrackIDMatchedEDepMap[edep_ptr->TrackID()].push_back(edep_ptr); + total_energy += edep_ptr->E(); + } + + // calculate EDep properties by TrackID + std::map TrackIDMatchedPropertyMap; + for (auto const& track_edeps : TrackIDMatchedEDepMap) + TrackIDMatchedPropertyMap[track_edeps.first] = CalcPropertiesFromEdeps(track_edeps.second, offset); + + // map it all out + std::map> EDepMatchedSubROIMap; // keys are indexes of edepPtrVec, values are vectors of indexes of subROIPropVec + std::map> TrackIDMatchedSubROIMap; // keys are TrackIDs, values are sets of indexes of subROIPropVec + std::map> SubROIMatchedEDepMap; // keys are indexes of subROIPropVec, values are vectors of indexes of edepPtrVec + std::map> SubROIMatchedTrackEnergyMap; // keys are indexes of subROIPropVec, values are maps of TrackIDs to matched energy (in MeV) + + // loop over EDeps + for (unsigned int i_e = 0; i_e < edepPtrVec.size(); ++i_e) + { + // get EDep properties + auto edep_ptr = edepPtrVec[i_e]; + const geo::TPCGeo& curTPCGeom = geometry->PositionToTPC(edep_ptr->MidPoint()); + double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now + double zeroTick = detPropData.ConvertXToTicks(0, curTPCGeom.Plane(0).ID()); + auto edep_tick = ticksPercm * edep_ptr->X() + (zeroTick + offset) + tickOffset; + edep_tick = detPropData.ConvertXToTicks(edep_ptr->X(), curTPCGeom.Plane(0).ID()) + offset + tickOffset; + + // loop over subROIs + unsigned int closest_hit = std::numeric_limits::max(); + float min_dist = std::numeric_limits::max(); + for (unsigned int i_h = 0; i_h < subROIPropVec.size(); ++i_h) + { + auto subroi_prop = subROIPropVec[i_h]; + if (edep_tick > subroi_prop.center-subroi_prop.sigma && edep_tick < subroi_prop.center+subroi_prop.sigma) + { + EDepMatchedSubROIMap[i_e].push_back(i_h); + TrackIDMatchedSubROIMap[edep_ptr->TrackID()].emplace(i_h); + } + float hit_dist = std::abs(edep_tick - subroi_prop.center) / subroi_prop.sigma; + if (hit_dist < min_dist) + { + closest_hit = i_h; + min_dist = hit_dist; + } + } + + // if EDep is less than 2.5 units away from closest subROI, assign it to that subROI + if (min_dist < 5) // try 5 for testing purposes + { + auto i_h = closest_hit; + SubROIMatchedEDepMap[i_h].push_back(i_e); + SubROIMatchedTrackEnergyMap[i_h][edep_ptr->TrackID()] += edep_ptr->E(); + } + } + + // convert to desired format (possibly a better way to do this...?) + std::map> ReturnMap; + for (auto it_h = SubROIMatchedEDepMap.begin(); it_h != SubROIMatchedEDepMap.end(); ++it_h) + { + auto key = subROIPropVec[it_h->first].key; + for (auto const& i_e : it_h->second) + { + ReturnMap[key].push_back(edepPtrVec[i_e]); + } + } + + return ReturnMap; +} + +//--- CalcPropertiesFromEdeps --- +sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEdeps(std::vector const& edepPtrVec, double offset) +{ + //split the edeps by TrackID + std::map< int, std::vector > edepptrs_by_trkid; + std::map< int, double > energy_per_trkid; + for(auto const& edep_ptr : edepPtrVec) + { + edepptrs_by_trkid[edep_ptr->TrackID()].push_back(edep_ptr); + energy_per_trkid[edep_ptr->TrackID()]+=edep_ptr->E(); + } + + int trkid_max = std::numeric_limits::min(); + double energy_max = std::numeric_limits::min(); + for(auto const& e_p_id : energy_per_trkid) + { + if(e_p_id.second > energy_max) + { + trkid_max = e_p_id.first; + energy_max = e_p_id.second; + } + } + + auto edepPtrVecMaxE = edepptrs_by_trkid[trkid_max]; + + //first, let's loop over all edeps and get an average weight scale... + sys::WireModUtility::TruthProperties_t edep_props; + double total_energy_all = 0.0; + sys::WireModUtility::ScaleValues_t scales_e_weighted[3]; + for(size_t i_p = 0; i_p < 3; ++i_p) + { + scales_e_weighted[i_p].r_Q = 0.0; + scales_e_weighted[i_p].r_sigma = 0.0; + } + + for(auto const edep_ptr : edepPtrVec) + { + if (edep_ptr->StepLength() == 0) + continue; + + edep_props.x = edep_ptr->X(); + edep_props.y = edep_ptr->Y(); + edep_props.z = edep_ptr->Z(); + + edep_props.dxdr = (edep_ptr->EndX() - edep_ptr->StartX()) / edep_ptr->StepLength(); + edep_props.dydr = (edep_ptr->EndY() - edep_ptr->StartY()) / edep_ptr->StepLength(); + edep_props.dzdr = (edep_ptr->EndZ() - edep_ptr->StartZ()) / edep_ptr->StepLength(); + edep_props.dedr = edep_ptr->E() / edep_ptr->StepLength(); + + total_energy_all += edep_ptr->E(); + + const geo::TPCGeo& curTPCGeom = geometry->PositionToTPC(edep_ptr->MidPoint()); + + for (size_t i_p = 0; i_p < 3; ++i_p) + { + auto scales = GetViewScaleValues(edep_props, curTPCGeom.Plane(i_p).View()); + scales_e_weighted[i_p].r_Q += edep_ptr->E()*scales.r_Q; + scales_e_weighted[i_p].r_sigma += edep_ptr->E()*scales.r_sigma; + } + } + + for(size_t i_p = 0; i_p < 3; ++i_p) + { + if (total_energy_all > 0) + { + scales_e_weighted[i_p].r_Q = scales_e_weighted[i_p].r_Q / total_energy_all; + scales_e_weighted[i_p].r_sigma = scales_e_weighted[i_p].r_sigma / total_energy_all; + } + if (scales_e_weighted[i_p].r_Q == 0) + scales_e_weighted[i_p].r_Q = 1; + if (scales_e_weighted[i_p].r_sigma == 0) + scales_e_weighted[i_p].r_sigma = 1; + } + + TruthProperties_t edep_col_properties; + + //copy in the scales that were calculated before + for(size_t i_p = 0; i_p < 3; ++i_p) + { + edep_col_properties.scales_avg[i_p].r_Q = scales_e_weighted[i_p].r_Q; + edep_col_properties.scales_avg[i_p].r_sigma = scales_e_weighted[i_p].r_sigma; + } + + // calculations happen here + edep_col_properties.x = 0.; + edep_col_properties.x_rms = 0.; + edep_col_properties.x_rms_noWeight = 0.; + edep_col_properties.x_min = std::numeric_limits::max(); + edep_col_properties.x_max = std::numeric_limits::min(); + + edep_col_properties.y = 0.; + edep_col_properties.z = 0.; + edep_col_properties.dxdr = 0.; + edep_col_properties.dydr = 0.; + edep_col_properties.dzdr = 0.; + edep_col_properties.dedr = 0.; + + double total_energy = 0.0; + for (auto const& edep_ptr : edepPtrVecMaxE) + { + edep_col_properties.x += edep_ptr->X()*edep_ptr->E(); + edep_col_properties.x_min = (edep_ptr->X() < edep_col_properties.x_min) ? edep_ptr->X() : edep_col_properties.x_min; + edep_col_properties.x_max = (edep_ptr->X() > edep_col_properties.x_max) ? edep_ptr->X() : edep_col_properties.x_max; + total_energy += edep_ptr->E(); + edep_col_properties.y += edep_ptr->Y()*edep_ptr->E(); + edep_col_properties.z += edep_ptr->Z()*edep_ptr->E(); + + if (edep_ptr->StepLength() == 0) + continue; + + edep_col_properties.dxdr += edep_ptr->E()*(edep_ptr->EndX() - edep_ptr->StartX()) / edep_ptr->StepLength(); + edep_col_properties.dydr += edep_ptr->E()*(edep_ptr->EndY() - edep_ptr->StartY()) / edep_ptr->StepLength(); + edep_col_properties.dzdr += edep_ptr->E()*(edep_ptr->EndZ() - edep_ptr->StartZ()) / edep_ptr->StepLength(); + edep_col_properties.dedr += edep_ptr->E()*edep_ptr->E() / edep_ptr->StepLength(); + } + + if (total_energy > 0) + { + edep_col_properties.x = edep_col_properties.x / total_energy; + edep_col_properties.y = edep_col_properties.y / total_energy; + edep_col_properties.z = edep_col_properties.z / total_energy; + edep_col_properties.dxdr = edep_col_properties.dxdr / total_energy; + edep_col_properties.dydr = edep_col_properties.dydr / total_energy; + edep_col_properties.dzdr = edep_col_properties.dzdr / total_energy; + edep_col_properties.dedr = edep_col_properties.dedr / total_energy; + } + + for (auto const& edep_ptr : edepPtrVecMaxE) + { + edep_col_properties.x_rms += (edep_ptr->X()-edep_col_properties.x)*(edep_ptr->X()-edep_col_properties.x)*edep_ptr->E(); + edep_col_properties.x_rms_noWeight += (edep_ptr->X()-edep_col_properties.x)*(edep_ptr->X()-edep_col_properties.x); + } + edep_col_properties.x_rms_noWeight = std::sqrt(edep_col_properties.x_rms_noWeight); + + if (total_energy > 0) + edep_col_properties.x_rms = std::sqrt(edep_col_properties.x_rms/total_energy); + + const geo::TPCGeo& tpcGeom = geometry->PositionToTPC({edep_col_properties.x, edep_col_properties.y, edep_col_properties.z}); + double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now + edep_col_properties.tick = detPropData.ConvertXToTicks(edep_col_properties.x , tpcGeom.Plane(0).ID()) + offset + tickOffset; + edep_col_properties.tick_rms = ticksPercm*edep_col_properties.x_rms; + edep_col_properties.tick_rms_noWeight = ticksPercm*edep_col_properties.x_rms_noWeight; + edep_col_properties.tick_min = detPropData.ConvertXToTicks(edep_col_properties.x_min, tpcGeom.Plane(0).ID()) + offset + tickOffset; + edep_col_properties.tick_max = detPropData.ConvertXToTicks(edep_col_properties.x_max, tpcGeom.Plane(0).ID()) + offset + tickOffset; + edep_col_properties.total_energy = total_energy; + + + return edep_col_properties; +} + +//--- GetScaleValues --- +sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetScaleValues(sys::WireModUtility::TruthProperties_t const& truth_props, sys::WireModUtility::ROIProperties_t const& roi_vals) +{ + sys::WireModUtility::ScaleValues_t scales; + sys::WireModUtility::ScaleValues_t channelScales = GetChannelScaleValues(truth_props, roi_vals.channel); + sys::WireModUtility::ScaleValues_t viewScales = GetViewScaleValues(truth_props, roi_vals.view); + scales.r_Q = channelScales.r_Q * viewScales.r_Q; + scales.r_sigma = channelScales.r_sigma * viewScales.r_sigma; + return scales; +} + +//--- GetChannelScaleValues --- +sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetChannelScaleValues(sys::WireModUtility::TruthProperties_t const& truth_props, raw::ChannelID_t const& channel) +{ + // initialize return + sys::WireModUtility::ScaleValues_t scales; + scales.r_Q = 1.0; + scales.r_sigma = 1.0; + + // try to get geo + // if not in a TPC return default values + double const truth_coords[3] = {truth_props.x, truth_props.y, truth_props.z}; + geo::TPCGeo const* curTPCGeomPtr = geometry->PositionToTPCptr(geo::vect::makePointFromCoords(truth_coords)); + if (curTPCGeomPtr == nullptr) + return scales; + + if (applyChannelScale) + { + if (spline_Charge_Channel == nullptr || + spline_Sigma_Channel == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply channel scale factor, but could not find splines. Check that you have set those in the utility."; + scales.r_Q *= spline_Charge_Channel->Eval(channel); + scales.r_sigma *= spline_Sigma_Channel ->Eval(channel); + } + + return scales; +} + +//--- GetViewScaleValues --- +sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys::WireModUtility::TruthProperties_t const& truth_props, geo::View_t const& view) +{ + // initialize return + sys::WireModUtility::ScaleValues_t scales; + scales.r_Q = 1.0; + scales.r_sigma = 1.0; + + double temp_scale=1.0; + + // try to get geo + // if not in a TPC return default values + double const truth_coords[3] = {truth_props.x, truth_props.y, truth_props.z}; + geo::TPCGeo const* curTPCGeomPtr = geometry->PositionToTPCptr(geo::vect::makePointFromCoords(truth_coords)); + if (curTPCGeomPtr == nullptr) + return scales; + + // get the plane number by the view + size_t plane = curTPCGeomPtr->Plane(view).ID().Plane; + + if (applyXScale) + { + if (splines_Charge_X[plane] == nullptr || + splines_Sigma_X [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply X scale factor, but could not find splines. Check that you have set those in the utility."; + scales.r_Q *= splines_Charge_X[plane]->Eval(truth_props.x); + scales.r_sigma *= splines_Sigma_X [plane]->Eval(truth_props.x); + } + + if (applyYZScale) + { + if (graph2Ds_Charge_YZ[plane] == nullptr || + 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); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = graph2Ds_Sigma_YZ [plane]->Interpolate(truth_props.z, truth_props.y); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + if (applyXZAngleScale) + { + if (splines_Charge_XZAngle[plane] == nullptr || + splines_Sigma_XZAngle [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply XZ-angle scale factor, but could not find splines. Check that you have set those in the utility."; + scales.r_Q *= splines_Charge_XZAngle[plane]->Eval(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, curTPCGeomPtr->Plane(plane).ThetaZ())); + scales.r_sigma *= splines_Sigma_XZAngle [plane]->Eval(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, curTPCGeomPtr->Plane(plane).ThetaZ())); + } + if (applyYZAngleScale) + { + if (splines_Charge_YZAngle[plane] == nullptr || + splines_Sigma_YZAngle [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply YZ-angle scale factor, but could not find splines. Check that you have set those in the utility."; + scales.r_Q *= splines_Charge_YZAngle[plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, curTPCGeomPtr->Plane(plane).ThetaZ())); + scales.r_sigma *= splines_Sigma_YZAngle [plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, curTPCGeomPtr->Plane(plane).ThetaZ())); + } + + if(applydEdXScale) + { + if (splines_Charge_dEdX[plane] == nullptr || + splines_Sigma_dEdX [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply dEdX scale factor, but could not find splines. Check that you have set those in the utility."; + scales.r_Q *= splines_Charge_dEdX[plane]->Eval(truth_props.dedr); + scales.r_sigma *= splines_Sigma_dEdX [plane]->Eval(truth_props.dedr); + } + + return scales; +} + +//--- ModifyROI --- +void sys::WireModUtility::ModifyROI(std::vector & roi_data, + sys::WireModUtility::ROIProperties_t const& roi_prop, + std::vector const& subROIPropVec, + std::map const& subROIScaleMap) +{ + // do you want a bunch of messages? + bool verbose = false; + + // initialize some values + double q_orig = 0.0; + double q_mod = 0.0; + double scale_ratio = 1.0; + + // loop over the ticks + for(size_t i_t = 0; i_t < roi_data.size(); ++i_t) + { + // reset your values + q_orig = 0.0; + q_mod = 0.0; + scale_ratio = 1.0; + + // 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); + + 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; + } + + // do some sanity checks + 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 { + 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 << "\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 + return; +} diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh new file mode 100644 index 000000000..30922d846 --- /dev/null +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -0,0 +1,218 @@ +#ifndef sbncode_WireMod_Utility_WireModUtility_hh +#define sbncode_WireMod_Utility_WireModUtility_hh + +//std includes +#include + +//ROOT includes +#include "TFile.h" +#include "TSpline.h" +#include "TGraph2D.h" +#include "TNtuple.h" + +//Framework includes +#include "larcorealg/Geometry/GeometryCore.h" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" +#include "lardataobj/RecoBase/Track.h" +#include "lardataobj/RecoBase/TrackHitMeta.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/Wire.h" +#include "lardataobj/Simulation/SimEnergyDeposit.h" +#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" +#include "nusimdata/SimulationBase/MCParticle.h" +#include "nusimdata/SimulationBase/MCTruth.h" + +#include +#include + +namespace sys { + class WireModUtility{ + // for now make everything public, though it's probably a good idea to think about what doesn't need to be + public: + // detector constants, should be set by geometry service + // the notes at the end refer to their old names in the MircoBooNE code which preceded this + // TODO: how best to initialize the splines/graphs? + const geo::GeometryCore* geometry; // save the TPC geometry + const detinfo::DetectorPropertiesData& detPropData; // save the detector property data + bool applyChannelScale; // do we scale with channel? + bool applyXScale; // do we scale with X? + bool applyYZScale; // do we scale with YZ? + bool applyXZAngleScale; // do we scale with XZ angle? + bool applyYZAngleScale; // do we scale with YZ angle? + bool applydEdXScale; // do we scale with dEdx? + double readoutWindowTicks; // how many ticks are in the readout window? + double tickOffset; // do we want an offset in the ticks? + + TSpline3* spline_Charge_Channel; // the spline for the charge correction by channel + TSpline3* spline_Sigma_Channel; // the spline for the width correction by channel + std::vector splines_Charge_X; // the splines for the charge correction in X + std::vector splines_Sigma_X; // the splines for the width correction in X + std::vector splines_Charge_XZAngle; // the splines for the charge correction in XZ angle + std::vector splines_Sigma_XZAngle; // the splines for the width correction in XZ angle + std::vector splines_Charge_YZAngle; // the splines for the charge correction in YZ angle + std::vector splines_Sigma_YZAngle; // the splines for the width correction in YZ angle + std::vector splines_Charge_dEdX; // the splines for the charge correction in dEdX + std::vector splines_Sigma_dEdX; // the splines for the width correction in dEdX + std::vector graph2Ds_Charge_YZ; // the graphs for the charge correction in YZ + std::vector graph2Ds_Sigma_YZ; // the graphs for the width correction in YZ + + // lets try making a constructor here + // assume we can get a geometry service, a detector clcok, and a detector properties + // pass the CryoStat and TPC IDs because it's IDs all the way down + // set some optional args fpr the booleans, the readout window, and the offset + WireModUtility(const geo::GeometryCore* geom, const detinfo::DetectorPropertiesData& detProp, + const bool& arg_ApplyChannelScale = false, + const bool& arg_ApplyXScale = true, + const bool& arg_ApplyYZScale = true, + const bool& arg_ApplyXZAngleScale = true, + const bool& arg_ApplyYZAngleScale = true, + const bool& arg_ApplydEdXScale = true, + const double& arg_TickOffset = 0) + : geometry(geom), + detPropData(detProp), + applyChannelScale(arg_ApplyChannelScale), + applyXScale(arg_ApplyXScale), + applyYZScale(arg_ApplyYZScale), + applyXZAngleScale(arg_ApplyXZAngleScale), + applyYZAngleScale(arg_ApplyYZAngleScale), + applydEdXScale(arg_ApplydEdXScale), + 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 + { + } + + // typedefs + typedef std::pair ROI_Key_t; + typedef std::pair SubROI_Key_t; + + typedef struct ROIProperties + { + ROI_Key_t key; + raw::ChannelID_t channel; + geo::View_t view; + float begin; + float end; + float total_q; + float center; //charge weighted center of ROI + float sigma; //charge weighted RMS of ROI + } ROIProperties_t; + + typedef struct SubROIProperties + { + SubROI_Key_t key; + raw::ChannelID_t channel; + geo::View_t view; + float total_q; + float center; + float sigma; + } SubROIProperties_t; + + typedef struct ScaleValues + { + double r_Q; + double r_sigma; + } ScaleValues_t; + + typedef struct TruthProperties + { + float x; + float x_rms; + float x_rms_noWeight; + float tick; + float tick_rms; + float tick_rms_noWeight; + float total_energy; + float x_min; + float x_max; + float tick_min; + float tick_max; + float y; + float z; + double dxdr; + double dydr; + double dzdr; + double dedr; + ScaleValues_t scales_avg[3]; + } TruthProperties_t; + + // internal containers + std::map< ROI_Key_t,std::vector > ROIMatchedEdepMap; + std::map< ROI_Key_t,std::vector > ROIMatchedHitMap; + + // some useful functions + // geometries + double planeXToTick(double xPos, int planeNumber, const geo::TPCGeo& tpcGeom, double offset = 0) { return detPropData.ConvertXToTicks(xPos, tpcGeom.Plane(planeNumber).ID()) + offset; } + + bool planeXInWindow(double xPos, int planeNumber, const geo::TPCGeo& tpcGeom, double offset = 0) + { + double tick = planeXToTick(xPos, planeNumber, tpcGeom, offset); + return (tick > 0 && tick <= detPropData.ReadOutWindowSize()); + } + + // for this function: in the future if we want to use non-gaussian functions make this take a vector of parameters + // the another wiremod utility could overwrite the ``fitFunc'' with some non-standard function + // would require a fair bit of remodling (ie q and sigma would need to be replace with, eg, funcVar[0] and funcVar[1] and probs a bunch of loops) + // so lets worry about that later + double gausFunc(double t, double mean, double sigma, double a = 1.0) + { + return (a / (sigma * std::sqrt(2 * util::pi()))) * std::exp(-0.5 * std::pow((t - mean)/sigma, 2)); + } + + double FoldAngle(double theta) + { + return (std::abs(theta) > 0.5 * util::pi()) ? util::pi() - std::abs(theta) : std::abs(theta); + } + + double ThetaXZ_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) + { + double planeAngleRad = planeAngle * (util::pi() / 180.0); + double sinPlaneAngle = std::sin(planeAngleRad); + double cosPlaneAngle = std::cos(planeAngleRad); + + //double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; // don't need to rotate Y for this angle + double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; + + double theta = std::atan2(dxdr, dzdrPlaneRel); + return FoldAngle(theta); + } + + double ThetaYZ_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) + { + double planeAngleRad = planeAngle * (util::pi() / 180.0); + double sinPlaneAngle = std::sin(planeAngleRad); + double cosPlaneAngle = std::cos(planeAngleRad); + + double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; + double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; + + double theta = std::atan2(dydrPlaneRel, dzdrPlaneRel); + return FoldAngle(theta); + } + + // theste are set in the .cc file + ROIProperties_t CalcROIProperties(recob::Wire const&, size_t const&); + + std::vector> GetTargetROIs(sim::SimEnergyDeposit const&, double offset); + std::vector> GetHitTargetROIs(recob::Hit const&); + + 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); + + TruthProperties_t CalcPropertiesFromEdeps(std::vector const&, double offset); + + 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&); + + void ModifyROI(std::vector &, + ROIProperties_t const &, + std::vector const&, + std::map const&); + }; // end class +} // end namespace + +#endif // sbncode_WireMod_Utility_WireModUtility_hh From f25dccc5c95881f1751f8a631e6346d016bfd6fe Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Wed, 19 Mar 2025 13:59:33 -0500 Subject: [PATCH 02/21] update for v10 --- sbncode/WireMod/Utility/WireModUtility.cc | 45 ++++++++++++++++------- sbncode/WireMod/Utility/WireModUtility.hh | 16 ++++++-- 2 files changed, 43 insertions(+), 18 deletions(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 71b0bb3ed..d61bac5ac 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -56,6 +56,20 @@ std::vector> sys::WireModUtility::GetTarge return target_roi_vec; // iterate over planes + for (auto const& plane : wireReadout->Iterate(curTPCGeomPtr->ID())) { + // check the wire exists in this plane + int wireNumber = int(0.5 + plane.WireCoordinate(shifted_edep.MidPoint())); + if ((wireNumber < 0) || (wireNumber >= (int)plane.Nwires())) + continue; + + // reconstruct the wireID from the position + geo::WireID edep_wireID = plane.NearestWireID(shifted_edep.MidPoint()); + + // if the deposition is inside the range where it would leave a signal on the wire, look for the ROI there + if (planeXInWindow(shifted_edep.X(), plane, *curTPCGeomPtr, offset + tickOffset)) + target_roi_vec.emplace_back(wireReadout->PlaneWireToChannel(edep_wireID), std::round(planeXToTick(shifted_edep.X(), plane, *curTPCGeomPtr, offset + tickOffset))); + } + /* for (size_t i_p = 0; i_p < curTPCGeomPtr->Nplanes(); ++i_p) { // check the wire exists in this plane @@ -70,6 +84,7 @@ std::vector> sys::WireModUtility::GetTarge if (planeXInWindow(shifted_edep.X(), i_p, *curTPCGeomPtr, offset + tickOffset)) target_roi_vec.emplace_back(geometry->PlaneWireToChannel(edep_wireID), std::round(planeXToTick(shifted_edep.X(), i_p, *curTPCGeomPtr, offset + tickOffset))); } + */ return target_roi_vec; } @@ -251,10 +266,11 @@ std::mapPositionToTPC(edep_ptr->MidPoint()); + const auto plane0 = wireReadout->FirstPlane(curTPCGeom.ID()); double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now - double zeroTick = detPropData.ConvertXToTicks(0, curTPCGeom.Plane(0).ID()); + double zeroTick = detPropData.ConvertXToTicks(0, plane0.ID()); auto edep_tick = ticksPercm * edep_ptr->X() + (zeroTick + offset) + tickOffset; - edep_tick = detPropData.ConvertXToTicks(edep_ptr->X(), curTPCGeom.Plane(0).ID()) + offset + tickOffset; + edep_tick = detPropData.ConvertXToTicks(edep_ptr->X(), plane0.ID()) + offset + tickOffset; // loop over subROIs unsigned int closest_hit = std::numeric_limits::max(); @@ -350,10 +366,9 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd total_energy_all += edep_ptr->E(); const geo::TPCGeo& curTPCGeom = geometry->PositionToTPC(edep_ptr->MidPoint()); - - for (size_t i_p = 0; i_p < 3; ++i_p) - { - auto scales = GetViewScaleValues(edep_props, curTPCGeom.Plane(i_p).View()); + for (auto const& plane : wireReadout->Iterate(curTPCGeom.ID())) { + int i_p = plane.ID().Plane; + auto scales = GetViewScaleValues(edep_props, plane.View()); scales_e_weighted[i_p].r_Q += edep_ptr->E()*scales.r_Q; scales_e_weighted[i_p].r_sigma += edep_ptr->E()*scales.r_sigma; } @@ -436,12 +451,13 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_col_properties.x_rms = std::sqrt(edep_col_properties.x_rms/total_energy); const geo::TPCGeo& tpcGeom = geometry->PositionToTPC({edep_col_properties.x, edep_col_properties.y, edep_col_properties.z}); + const auto plane0 = wireReadout->FirstPlane(tpcGeom.ID()); double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now - edep_col_properties.tick = detPropData.ConvertXToTicks(edep_col_properties.x , tpcGeom.Plane(0).ID()) + offset + tickOffset; + edep_col_properties.tick = detPropData.ConvertXToTicks(edep_col_properties.x , plane0.ID()) + offset + tickOffset; edep_col_properties.tick_rms = ticksPercm*edep_col_properties.x_rms; edep_col_properties.tick_rms_noWeight = ticksPercm*edep_col_properties.x_rms_noWeight; - edep_col_properties.tick_min = detPropData.ConvertXToTicks(edep_col_properties.x_min, tpcGeom.Plane(0).ID()) + offset + tickOffset; - edep_col_properties.tick_max = detPropData.ConvertXToTicks(edep_col_properties.x_max, tpcGeom.Plane(0).ID()) + offset + tickOffset; + edep_col_properties.tick_min = detPropData.ConvertXToTicks(edep_col_properties.x_min, plane0.ID()) + offset + tickOffset; + edep_col_properties.tick_max = detPropData.ConvertXToTicks(edep_col_properties.x_max, plane0.ID()) + offset + tickOffset; edep_col_properties.total_energy = total_energy; @@ -505,7 +521,8 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: return scales; // get the plane number by the view - size_t plane = curTPCGeomPtr->Plane(view).ID().Plane; + auto const& plane_obj = wireReadout->Plane(curTPCGeomPtr->ID(), view); + unsigned int plane = plane_obj.ID().Plane; if (applyXScale) { @@ -536,8 +553,8 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: splines_Sigma_XZAngle [plane] == nullptr ) throw cet::exception("WireModUtility") << "Tried to apply XZ-angle scale factor, but could not find splines. Check that you have set those in the utility."; - scales.r_Q *= splines_Charge_XZAngle[plane]->Eval(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, curTPCGeomPtr->Plane(plane).ThetaZ())); - scales.r_sigma *= splines_Sigma_XZAngle [plane]->Eval(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, curTPCGeomPtr->Plane(plane).ThetaZ())); + scales.r_Q *= splines_Charge_XZAngle[plane]->Eval(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + scales.r_sigma *= splines_Sigma_XZAngle [plane]->Eval(ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); } if (applyYZAngleScale) { @@ -545,8 +562,8 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: splines_Sigma_YZAngle [plane] == nullptr ) throw cet::exception("WireModUtility") << "Tried to apply YZ-angle scale factor, but could not find splines. Check that you have set those in the utility."; - scales.r_Q *= splines_Charge_YZAngle[plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, curTPCGeomPtr->Plane(plane).ThetaZ())); - scales.r_sigma *= splines_Sigma_YZAngle [plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, curTPCGeomPtr->Plane(plane).ThetaZ())); + scales.r_Q *= splines_Charge_YZAngle[plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + scales.r_sigma *= splines_Sigma_YZAngle [plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); } if(applydEdXScale) diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index 30922d846..570d262ec 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -12,6 +12,7 @@ //Framework includes #include "larcorealg/Geometry/GeometryCore.h" +#include "larcorealg/Geometry/WireReadoutGeom.h" #include "lardataalg/DetectorInfo/DetectorPropertiesData.h" #include "lardataobj/RecoBase/Track.h" #include "lardataobj/RecoBase/TrackHitMeta.h" @@ -33,6 +34,7 @@ namespace sys { // the notes at the end refer to their old names in the MircoBooNE code which preceded this // TODO: how best to initialize the splines/graphs? const geo::GeometryCore* geometry; // save the TPC geometry + const geo::WireReadoutGeom* wireReadout; // new for LarSoft v10 const detinfo::DetectorPropertiesData& detPropData; // save the detector property data bool applyChannelScale; // do we scale with channel? bool applyXScale; // do we scale with X? @@ -60,7 +62,9 @@ namespace sys { // assume we can get a geometry service, a detector clcok, and a detector properties // pass the CryoStat and TPC IDs because it's IDs all the way down // set some optional args fpr the booleans, the readout window, and the offset - WireModUtility(const geo::GeometryCore* geom, const detinfo::DetectorPropertiesData& detProp, + WireModUtility(const geo::GeometryCore* geom, + const geo::WireReadoutGeom* wireRead, + const detinfo::DetectorPropertiesData& detProp, const bool& arg_ApplyChannelScale = false, const bool& arg_ApplyXScale = true, const bool& arg_ApplyYZScale = true, @@ -69,6 +73,7 @@ namespace sys { const bool& arg_ApplydEdXScale = true, const double& arg_TickOffset = 0) : geometry(geom), + wireReadout(wireRead), detPropData(detProp), applyChannelScale(arg_ApplyChannelScale), applyXScale(arg_ApplyXScale), @@ -141,11 +146,14 @@ namespace sys { // some useful functions // geometries - double planeXToTick(double xPos, int planeNumber, const geo::TPCGeo& tpcGeom, double offset = 0) { return detPropData.ConvertXToTicks(xPos, tpcGeom.Plane(planeNumber).ID()) + offset; } + // TODO is this the most efficient for new v10 iterators? + double planeXToTick(double xPos, const geo::PlaneGeo& plane, const geo::TPCGeo& tpcGeom, double offset = 0) { + return detPropData.ConvertXToTicks(xPos, plane.ID()) + offset; + } - bool planeXInWindow(double xPos, int planeNumber, const geo::TPCGeo& tpcGeom, double offset = 0) + bool planeXInWindow(double xPos, const geo::PlaneGeo& plane, const geo::TPCGeo& tpcGeom, double offset = 0) { - double tick = planeXToTick(xPos, planeNumber, tpcGeom, offset); + double tick = planeXToTick(xPos, plane, tpcGeom, offset); return (tick > 0 && tick <= detPropData.ReadOutWindowSize()); } From 52892ebef281402d0068de0d935b7137292459c0 Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Thu, 16 Oct 2025 17:36:38 -0500 Subject: [PATCH 03/21] WireMod for sbnd production branch --- sbncode/WireMod/Utility/WireModUtility.cc | 45 ++++++++++++++++++++++- sbncode/WireMod/Utility/WireModUtility.hh | 32 ++++++++++++++++ 2 files changed, 76 insertions(+), 1 deletion(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index d61bac5ac..126044753 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -361,6 +361,7 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_props.dxdr = (edep_ptr->EndX() - edep_ptr->StartX()) / edep_ptr->StepLength(); edep_props.dydr = (edep_ptr->EndY() - edep_ptr->StartY()) / edep_ptr->StepLength(); edep_props.dzdr = (edep_ptr->EndZ() - edep_ptr->StartZ()) / edep_ptr->StepLength(); + edep_props.dedr = edep_ptr->E() / edep_ptr->StepLength(); total_energy_all += edep_ptr->E(); @@ -408,6 +409,7 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_col_properties.dxdr = 0.; edep_col_properties.dydr = 0.; edep_col_properties.dzdr = 0.; + edep_col_properties.dedr = 0.; double total_energy = 0.0; @@ -426,6 +428,7 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_col_properties.dxdr += edep_ptr->E()*(edep_ptr->EndX() - edep_ptr->StartX()) / edep_ptr->StepLength(); edep_col_properties.dydr += edep_ptr->E()*(edep_ptr->EndY() - edep_ptr->StartY()) / edep_ptr->StepLength(); edep_col_properties.dzdr += edep_ptr->E()*(edep_ptr->EndZ() - edep_ptr->StartZ()) / edep_ptr->StepLength(); + edep_col_properties.dedr += edep_ptr->E()*edep_ptr->E() / edep_ptr->StepLength(); } @@ -437,6 +440,7 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_col_properties.dxdr = edep_col_properties.dxdr / total_energy; edep_col_properties.dydr = edep_col_properties.dydr / total_energy; edep_col_properties.dzdr = edep_col_properties.dzdr / total_energy; + edep_col_properties.dedr = edep_col_properties.dedr / total_energy; } @@ -566,7 +570,7 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: scales.r_sigma *= splines_Sigma_YZAngle [plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); } - if(applydEdXScale) + if (applydEdXScale) { if (splines_Charge_dEdX[plane] == nullptr || splines_Sigma_dEdX [plane] == nullptr ) @@ -576,6 +580,45 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: scales.r_sigma *= splines_Sigma_dEdX [plane]->Eval(truth_props.dedr); } + if (applyXXZAngleScale) + { + if (graph2Ds_Charge_XXZAngle[plane] == nullptr || + 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())); + 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())); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + if (applyXdQdXScale) + { + if (graph2Ds_Charge_XdQdX[plane] == nullptr || + 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); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = graph2Ds_Sigma_XdQdX [plane]->Interpolate(truth_props.x, truth_props.dqdr); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + if (applyXZAngledQdXScale) + { + if (graph2Ds_Charge_XZAngledQdX[plane] == nullptr || + 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); + 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); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + return scales; } diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index 570d262ec..6d7c1e31a 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -42,6 +42,12 @@ namespace sys { bool applyXZAngleScale; // do we scale with XZ angle? bool applyYZAngleScale; // do we scale with YZ angle? bool applydEdXScale; // do we scale with dEdx? +<<<<<<< HEAD +======= + bool applyXXZAngleScale; // do we scale with X vs XZ angle? + bool applyXdQdXScale; // do we scale with X vs dQ/dX? + bool applyXZAngledQdXScale; // do we scale with XZ angle vs dQ/dX? +>>>>>>> feature/hhausner_wiremod_v10_MoreScalingEnabled double readoutWindowTicks; // how many ticks are in the readout window? double tickOffset; // do we want an offset in the ticks? @@ -58,6 +64,16 @@ namespace sys { std::vector graph2Ds_Charge_YZ; // the graphs for the charge correction in YZ std::vector graph2Ds_Sigma_YZ; // the graphs for the width correction in YZ +<<<<<<< HEAD +======= + std::vector graph2Ds_Charge_XXZAngle; // the graphs for the charge correction in X vs XZ angle + std::vector graph2Ds_Sigma_XXZAngle; // the graphs for the width correction in X vs XZ angle + std::vector graph2Ds_Charge_XdQdX; // the graphs for charge correction in X vs dQ/dX + std::vector graph2Ds_Sigma_XdQdX; // the graphs for width correction in X vs dQ/dX + std::vector graph2Ds_Charge_XZAngledQdX; // the graphs for charge correction in XZ angle vs dQ/dX + std::vector graph2Ds_Sigma_XZAngledQdX; // the graphs for width correction in XZ angle vs dQ/dX + +>>>>>>> feature/hhausner_wiremod_v10_MoreScalingEnabled // lets try making a constructor here // assume we can get a geometry service, a detector clcok, and a detector properties // pass the CryoStat and TPC IDs because it's IDs all the way down @@ -71,6 +87,12 @@ namespace sys { const bool& arg_ApplyXZAngleScale = true, const bool& arg_ApplyYZAngleScale = true, const bool& arg_ApplydEdXScale = true, +<<<<<<< HEAD +======= + const bool& arg_ApplyXXZAngleScale = false, + const bool& arg_ApplyXdQdXScale = false, + const bool& arg_ApplyXZAngledQdXScale = false, +>>>>>>> feature/hhausner_wiremod_v10_MoreScalingEnabled const double& arg_TickOffset = 0) : geometry(geom), wireReadout(wireRead), @@ -81,6 +103,12 @@ namespace sys { applyXZAngleScale(arg_ApplyXZAngleScale), applyYZAngleScale(arg_ApplyYZAngleScale), applydEdXScale(arg_ApplydEdXScale), +<<<<<<< HEAD +======= + applyXXZAngleScale(arg_ApplyXXZAngleScale), + applyXdQdXScale(arg_ApplyXdQdXScale), + applyXZAngledQdXScale(arg_ApplyXZAngledQdXScale), +>>>>>>> feature/hhausner_wiremod_v10_MoreScalingEnabled 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 { @@ -136,6 +164,10 @@ namespace sys { double dxdr; double dydr; double dzdr; +<<<<<<< HEAD +======= + double dqdr; +>>>>>>> feature/hhausner_wiremod_v10_MoreScalingEnabled double dedr; ScaleValues_t scales_avg[3]; } TruthProperties_t; From cb5d01be0c19bee20c5b142bcfe2a8c9c3b51bfd Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Thu, 23 Oct 2025 12:01:41 -0500 Subject: [PATCH 04/21] fix build --- sbncode/WireMod/Utility/WireModUtility.hh | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index 6d7c1e31a..1ae74c152 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -42,12 +42,9 @@ namespace sys { bool applyXZAngleScale; // do we scale with XZ angle? bool applyYZAngleScale; // do we scale with YZ angle? bool applydEdXScale; // do we scale with dEdx? -<<<<<<< HEAD -======= bool applyXXZAngleScale; // do we scale with X vs XZ angle? bool applyXdQdXScale; // do we scale with X vs dQ/dX? bool applyXZAngledQdXScale; // do we scale with XZ angle vs dQ/dX? ->>>>>>> feature/hhausner_wiremod_v10_MoreScalingEnabled double readoutWindowTicks; // how many ticks are in the readout window? double tickOffset; // do we want an offset in the ticks? @@ -64,8 +61,6 @@ namespace sys { std::vector graph2Ds_Charge_YZ; // the graphs for the charge correction in YZ std::vector graph2Ds_Sigma_YZ; // the graphs for the width correction in YZ -<<<<<<< HEAD -======= std::vector graph2Ds_Charge_XXZAngle; // the graphs for the charge correction in X vs XZ angle std::vector graph2Ds_Sigma_XXZAngle; // the graphs for the width correction in X vs XZ angle std::vector graph2Ds_Charge_XdQdX; // the graphs for charge correction in X vs dQ/dX @@ -73,7 +68,6 @@ namespace sys { std::vector graph2Ds_Charge_XZAngledQdX; // the graphs for charge correction in XZ angle vs dQ/dX std::vector graph2Ds_Sigma_XZAngledQdX; // the graphs for width correction in XZ angle vs dQ/dX ->>>>>>> feature/hhausner_wiremod_v10_MoreScalingEnabled // lets try making a constructor here // assume we can get a geometry service, a detector clcok, and a detector properties // pass the CryoStat and TPC IDs because it's IDs all the way down @@ -87,12 +81,9 @@ namespace sys { const bool& arg_ApplyXZAngleScale = true, const bool& arg_ApplyYZAngleScale = true, const bool& arg_ApplydEdXScale = true, -<<<<<<< HEAD -======= const bool& arg_ApplyXXZAngleScale = false, const bool& arg_ApplyXdQdXScale = false, const bool& arg_ApplyXZAngledQdXScale = false, ->>>>>>> feature/hhausner_wiremod_v10_MoreScalingEnabled const double& arg_TickOffset = 0) : geometry(geom), wireReadout(wireRead), @@ -103,12 +94,9 @@ namespace sys { applyXZAngleScale(arg_ApplyXZAngleScale), applyYZAngleScale(arg_ApplyYZAngleScale), applydEdXScale(arg_ApplydEdXScale), -<<<<<<< HEAD -======= applyXXZAngleScale(arg_ApplyXXZAngleScale), applyXdQdXScale(arg_ApplyXdQdXScale), applyXZAngledQdXScale(arg_ApplyXZAngledQdXScale), ->>>>>>> feature/hhausner_wiremod_v10_MoreScalingEnabled 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 { @@ -164,10 +152,7 @@ namespace sys { double dxdr; double dydr; double dzdr; -<<<<<<< HEAD -======= double dqdr; ->>>>>>> feature/hhausner_wiremod_v10_MoreScalingEnabled double dedr; ScaleValues_t scales_avg[3]; } TruthProperties_t; From f8a9ea7209f516902d04e9c2c43ac0074fa270b8 Mon Sep 17 00:00:00 2001 From: hausnerh Date: Tue, 13 Jan 2026 10:39:54 -0600 Subject: [PATCH 05/21] add WireMod in X vs theta_XW --- sbncode/CAFMaker/RecoUtils/CMakeLists.txt | 1 + sbncode/WireMod/Utility/WireModUtility.cc | 63 +++++++---------------- sbncode/WireMod/Utility/WireModUtility.hh | 53 +++++++++++++------ 3 files changed, 56 insertions(+), 61 deletions(-) diff --git a/sbncode/CAFMaker/RecoUtils/CMakeLists.txt b/sbncode/CAFMaker/RecoUtils/CMakeLists.txt index b6abdbdb9..5401c9bbd 100644 --- a/sbncode/CAFMaker/RecoUtils/CMakeLists.txt +++ b/sbncode/CAFMaker/RecoUtils/CMakeLists.txt @@ -2,6 +2,7 @@ art_make_library( LIBRARY_NAME caf_RecoUtils SOURCE RecoUtils.cc LIBRARIES + sbnanaobj::StandardRecord art::Framework_Core art::Framework_Services_Registry art::Framework_Principal diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 126044753..0e4dc4da4 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -108,6 +108,11 @@ std::vector> sys::WireModUtility::GetHitTa return target_roi_vec; } +//--- FillROIMatchedIDEMap --- +void sys::WireModUtility::FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, double offset) +{ +} + //--- FillROIMatchedEdepMap --- void sys::WireModUtility::FillROIMatchedEdepMap(std::vector const& edepVec, std::vector const& wireVec, double offset) { @@ -361,7 +366,6 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_props.dxdr = (edep_ptr->EndX() - edep_ptr->StartX()) / edep_ptr->StepLength(); edep_props.dydr = (edep_ptr->EndY() - edep_ptr->StartY()) / edep_ptr->StepLength(); edep_props.dzdr = (edep_ptr->EndZ() - edep_ptr->StartZ()) / edep_ptr->StepLength(); - edep_props.dedr = edep_ptr->E() / edep_ptr->StepLength(); total_energy_all += edep_ptr->E(); @@ -409,7 +413,6 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_col_properties.dxdr = 0.; edep_col_properties.dydr = 0.; edep_col_properties.dzdr = 0.; - edep_col_properties.dedr = 0.; double total_energy = 0.0; @@ -428,7 +431,6 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_col_properties.dxdr += edep_ptr->E()*(edep_ptr->EndX() - edep_ptr->StartX()) / edep_ptr->StepLength(); edep_col_properties.dydr += edep_ptr->E()*(edep_ptr->EndY() - edep_ptr->StartY()) / edep_ptr->StepLength(); edep_col_properties.dzdr += edep_ptr->E()*(edep_ptr->EndZ() - edep_ptr->StartZ()) / edep_ptr->StepLength(); - edep_col_properties.dedr += edep_ptr->E()*edep_ptr->E() / edep_ptr->StepLength(); } @@ -440,7 +442,6 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_col_properties.dxdr = edep_col_properties.dxdr / total_energy; edep_col_properties.dydr = edep_col_properties.dydr / total_energy; edep_col_properties.dzdr = edep_col_properties.dzdr / total_energy; - edep_col_properties.dedr = edep_col_properties.dedr / total_energy; } @@ -551,6 +552,19 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: if(temp_scale>0.001) scales.r_sigma *= temp_scale; } + if (applyXXWScale) + { + if (graph2Ds_Charge_XXW[plane] == nullptr || + graph2Ds_Sigma_XXW [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply XXW scale factor, but could not find graphs. Check that you have set those in the utility."; + temp_scale = graph2Ds_Charge_XXW[plane]->Interpolate(truth_props.x, ThetaXW(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = graph2Ds_Sigma_XXW [plane]->Interpolate(truth_props.x, ThetaXW(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + if (applyXZAngleScale) { if (splines_Charge_XZAngle[plane] == nullptr || @@ -570,7 +584,7 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: scales.r_sigma *= splines_Sigma_YZAngle [plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); } - if (applydEdXScale) + if(applydEdXScale) { if (splines_Charge_dEdX[plane] == nullptr || splines_Sigma_dEdX [plane] == nullptr ) @@ -580,45 +594,6 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: scales.r_sigma *= splines_Sigma_dEdX [plane]->Eval(truth_props.dedr); } - if (applyXXZAngleScale) - { - if (graph2Ds_Charge_XXZAngle[plane] == nullptr || - 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())); - 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())); - if(temp_scale>0.001) scales.r_sigma *= temp_scale; - } - - if (applyXdQdXScale) - { - if (graph2Ds_Charge_XdQdX[plane] == nullptr || - 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); - if(temp_scale>0.001) scales.r_Q *= temp_scale; - - temp_scale = graph2Ds_Sigma_XdQdX [plane]->Interpolate(truth_props.x, truth_props.dqdr); - if(temp_scale>0.001) scales.r_sigma *= temp_scale; - } - - if (applyXZAngledQdXScale) - { - if (graph2Ds_Charge_XZAngledQdX[plane] == nullptr || - 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); - 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); - if(temp_scale>0.001) scales.r_sigma *= temp_scale; - } - return scales; } diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index 1ae74c152..0bc5c2d75 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -19,6 +19,7 @@ #include "lardataobj/RecoBase/Hit.h" #include "lardataobj/RecoBase/Wire.h" #include "lardataobj/Simulation/SimEnergyDeposit.h" +#include "lardataobj/Simulation/SimChannel.h" #include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" #include "nusimdata/SimulationBase/MCParticle.h" #include "nusimdata/SimulationBase/MCTruth.h" @@ -42,9 +43,7 @@ namespace sys { bool applyXZAngleScale; // do we scale with XZ angle? bool applyYZAngleScale; // do we scale with YZ angle? bool applydEdXScale; // do we scale with dEdx? - bool applyXXZAngleScale; // do we scale with X vs XZ angle? - bool applyXdQdXScale; // do we scale with X vs dQ/dX? - bool applyXZAngledQdXScale; // do we scale with XZ angle vs dQ/dX? + bool applyXXWScale; // do we scale with X vs ThXW? double readoutWindowTicks; // how many ticks are in the readout window? double tickOffset; // do we want an offset in the ticks? @@ -60,13 +59,8 @@ namespace sys { std::vector splines_Sigma_dEdX; // the splines for the width correction in dEdX std::vector graph2Ds_Charge_YZ; // the graphs for the charge correction in YZ std::vector graph2Ds_Sigma_YZ; // the graphs for the width correction in YZ - - std::vector graph2Ds_Charge_XXZAngle; // the graphs for the charge correction in X vs XZ angle - std::vector graph2Ds_Sigma_XXZAngle; // the graphs for the width correction in X vs XZ angle - std::vector graph2Ds_Charge_XdQdX; // the graphs for charge correction in X vs dQ/dX - std::vector graph2Ds_Sigma_XdQdX; // the graphs for width correction in X vs dQ/dX - std::vector graph2Ds_Charge_XZAngledQdX; // the graphs for charge correction in XZ angle vs dQ/dX - std::vector graph2Ds_Sigma_XZAngledQdX; // the graphs for width correction in XZ angle vs dQ/dX + std::vector graph2Ds_Charge_XXW; // the graphs for the charge correction in XXW + std::vector graph2Ds_Sigma_XXW; // the graphs for the width correction in XXW // lets try making a constructor here // assume we can get a geometry service, a detector clcok, and a detector properties @@ -81,9 +75,7 @@ namespace sys { const bool& arg_ApplyXZAngleScale = true, const bool& arg_ApplyYZAngleScale = true, const bool& arg_ApplydEdXScale = true, - const bool& arg_ApplyXXZAngleScale = false, - const bool& arg_ApplyXdQdXScale = false, - const bool& arg_ApplyXZAngledQdXScale = false, + const bool& arg_ApplyXXWScale = true, const double& arg_TickOffset = 0) : geometry(geom), wireReadout(wireRead), @@ -94,9 +86,7 @@ namespace sys { applyXZAngleScale(arg_ApplyXZAngleScale), applyYZAngleScale(arg_ApplyYZAngleScale), applydEdXScale(arg_ApplydEdXScale), - applyXXZAngleScale(arg_ApplyXXZAngleScale), - applyXdQdXScale(arg_ApplyXdQdXScale), - applyXZAngledQdXScale(arg_ApplyXZAngledQdXScale), + applyXXWScale(arg_ApplyXXWScale), 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 { @@ -152,7 +142,6 @@ namespace sys { double dxdr; double dydr; double dzdr; - double dqdr; double dedr; ScaleValues_t scales_avg[3]; } TruthProperties_t; @@ -213,6 +202,35 @@ namespace sys { double theta = std::atan2(dydrPlaneRel, dzdrPlaneRel); return FoldAngle(theta); } + + double ThetaXY_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) + { + double planeAngleRad = planeAngle * (util::pi() / 180.0); + double sinPlaneAngle = std::sin(planeAngleRad); + double cosPlaneAngle = std::cos(planeAngleRad); + + double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; + //double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; // don't need to rotate Z for this angle + + double theta = std::atan2(dxdr, dydrPlaneRel); + return FoldAngle(theta); + } + + double ThetaXW(double dxdr, double dydr, double dzdr, double planeAngle) + { + double planeAngleRad = planeAngle * (util::pi() / 180.0); + double sinPlaneAngle = std::sin(planeAngleRad); + double cosPlaneAngle = std::cos(planeAngleRad); + + double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; + //double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; // don't need to rotate Z for this angle + + //double theta = std::atan2(dxdr, dydrPlaneRel); + //return theta; + + double theta = std::atan(dydrPlaneRel / dxdr); + return FoldAngle(theta); + } // theste are set in the .cc file ROIProperties_t CalcROIProperties(recob::Wire const&, size_t const&); @@ -222,6 +240,7 @@ 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); std::vector CalcSubROIProperties(ROIProperties_t const&, std::vector const&); From eda854247bf749910de4716c181eaffb28b2ccad Mon Sep 17 00:00:00 2001 From: hausnerh Date: Tue, 3 Feb 2026 09:30:43 -0600 Subject: [PATCH 06/21] I believe this is the correct ThXW calculation... --- sbncode/WireMod/Utility/WireModUtility.hh | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index 0bc5c2d75..c51f35e15 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -222,14 +222,9 @@ namespace sys { double sinPlaneAngle = std::sin(planeAngleRad); double cosPlaneAngle = std::cos(planeAngleRad); - double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; - //double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; // don't need to rotate Z for this angle - - //double theta = std::atan2(dxdr, dydrPlaneRel); - //return theta; - - double theta = std::atan(dydrPlaneRel / dxdr); - return FoldAngle(theta); + double cosG = std::abs(dydr * sinPlaneAngle + dzdr * cosPlaneAngle); + double theta = std::atan(dxdr / cosG); + return std::abs(theta); } // theste are set in the .cc file From 60d87e6da249bf8d6ebdadb2391c6d209b9bdf9b Mon Sep 17 00:00:00 2001 From: hausnerh Date: Mon, 9 Mar 2026 14:22:10 -0500 Subject: [PATCH 07/21] Slight tweak to ThXW so it's calc'd once --- sbncode/WireMod/Utility/WireModUtility.cc | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 0e4dc4da4..0d56a24fa 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -558,10 +558,11 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: graph2Ds_Sigma_XXW [plane] == nullptr ) throw cet::exception("WireModUtility") << "Tried to apply XXW scale factor, but could not find graphs. Check that you have set those in the utility."; - temp_scale = graph2Ds_Charge_XXW[plane]->Interpolate(truth_props.x, ThetaXW(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + 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); if(temp_scale>0.001) scales.r_Q *= temp_scale; - temp_scale = graph2Ds_Sigma_XXW [plane]->Interpolate(truth_props.x, ThetaXW(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); + temp_scale = graph2Ds_Sigma_XXW [plane]->Interpolate(truth_props.x, thXW); if(temp_scale>0.001) scales.r_sigma *= temp_scale; } From 2060d5d4c738c7ea988df0dc3cecdb5df81345f5 Mon Sep 17 00:00:00 2001 From: hausnerh Date: Thu, 26 Mar 2026 11:11:49 -0500 Subject: [PATCH 08/21] bug fixes to the WM util + slightly smarting handling to the EDep projection --- sbncode/WireMod/Utility/WireModUtility.cc | 25 ++++++++++++++--------- sbncode/WireMod/Utility/WireModUtility.hh | 24 +++++++++------------- 2 files changed, 25 insertions(+), 24 deletions(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 0d56a24fa..3bd0834ed 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -254,11 +254,6 @@ std::mapE(); } - // calculate EDep properties by TrackID - std::map TrackIDMatchedPropertyMap; - for (auto const& track_edeps : TrackIDMatchedEDepMap) - TrackIDMatchedPropertyMap[track_edeps.first] = CalcPropertiesFromEdeps(track_edeps.second, offset); - // map it all out std::map> EDepMatchedSubROIMap; // keys are indexes of edepPtrVec, values are vectors of indexes of subROIPropVec std::map> TrackIDMatchedSubROIMap; // keys are TrackIDs, values are sets of indexes of subROIPropVec @@ -271,11 +266,18 @@ std::mapPositionToTPC(edep_ptr->MidPoint()); - const auto plane0 = wireReadout->FirstPlane(curTPCGeom.ID()); - double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now - double zeroTick = detPropData.ConvertXToTicks(0, plane0.ID()); - auto edep_tick = ticksPercm * edep_ptr->X() + (zeroTick + offset) + tickOffset; - edep_tick = detPropData.ConvertXToTicks(edep_ptr->X(), plane0.ID()) + offset + tickOffset; + std::map planeProjEdepTick; + for (const auto& planeGeom : wireReadout->Iterate(curTPCGeom.ID())) + { + double projTick = detPropData.ConvertXToTicks(edep_ptr->X(), planeGeom.ID()) + offset + tickOffset; + readout::ROPID ropID = wireReadout->WirePlaneToROP(planeGeom.ID()); + if (planeProjEdepTick.count(ropID) == 1) + { + std::cout << "MatchEdepsToSubROIs - Warning, overwriting tick projection at" << ropID.toString() + << "\n Changing from " << planeProjEdepTick[ropID] << " to " << projTick << std::endl; + } + planeProjEdepTick[ropID] = projTick; + } // loop over subROIs unsigned int closest_hit = std::numeric_limits::max(); @@ -283,6 +285,8 @@ std::mapChannelToROP(subroi_prop.channel); + double edep_tick = planeProjEdepTick.at(subroiROP); if (edep_tick > subroi_prop.center-subroi_prop.sigma && edep_tick < subroi_prop.center+subroi_prop.sigma) { EDepMatchedSubROIMap[i_e].push_back(i_h); @@ -455,6 +459,7 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd if (total_energy > 0) edep_col_properties.x_rms = std::sqrt(edep_col_properties.x_rms/total_energy); + // projecting edep to plane0, so be mindful when accessing the tick value const geo::TPCGeo& tpcGeom = geometry->PositionToTPC({edep_col_properties.x, edep_col_properties.y, edep_col_properties.z}); const auto plane0 = wireReadout->FirstPlane(tpcGeom.ID()); double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index c51f35e15..113bd9122 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -135,8 +135,8 @@ namespace sys { float total_energy; float x_min; float x_max; - float tick_min; - float tick_max; + float tick_min; // On Plane0!! + float tick_max; // On Plane0!! float y; float z; double dxdr; @@ -179,9 +179,8 @@ namespace sys { double ThetaXZ_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) { - double planeAngleRad = planeAngle * (util::pi() / 180.0); - double sinPlaneAngle = std::sin(planeAngleRad); - double cosPlaneAngle = std::cos(planeAngleRad); + double sinPlaneAngle = std::sin(planeAngle); + double cosPlaneAngle = std::cos(planeAngle); //double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; // don't need to rotate Y for this angle double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; @@ -192,9 +191,8 @@ namespace sys { double ThetaYZ_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) { - double planeAngleRad = planeAngle * (util::pi() / 180.0); - double sinPlaneAngle = std::sin(planeAngleRad); - double cosPlaneAngle = std::cos(planeAngleRad); + double sinPlaneAngle = std::sin(planeAngle); + double cosPlaneAngle = std::cos(planeAngle); double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; @@ -205,9 +203,8 @@ namespace sys { double ThetaXY_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) { - double planeAngleRad = planeAngle * (util::pi() / 180.0); - double sinPlaneAngle = std::sin(planeAngleRad); - double cosPlaneAngle = std::cos(planeAngleRad); + double sinPlaneAngle = std::sin(planeAngle); + double cosPlaneAngle = std::cos(planeAngle); double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; //double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; // don't need to rotate Z for this angle @@ -218,9 +215,8 @@ namespace sys { double ThetaXW(double dxdr, double dydr, double dzdr, double planeAngle) { - double planeAngleRad = planeAngle * (util::pi() / 180.0); - double sinPlaneAngle = std::sin(planeAngleRad); - double cosPlaneAngle = std::cos(planeAngleRad); + double sinPlaneAngle = std::sin(planeAngle); + double cosPlaneAngle = std::cos(planeAngle); double cosG = std::abs(dydr * sinPlaneAngle + dzdr * cosPlaneAngle); double theta = std::atan(dxdr / cosG); From ed07bb3d99bd0867273dfd74c4b29ac4fa44d807 Mon Sep 17 00:00:00 2001 From: hausnerh Date: Mon, 30 Mar 2026 15:13:21 -0500 Subject: [PATCH 09/21] fix a few spaces & typos --- sbncode/WireMod/Utility/WireModUtility.hh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index 113bd9122..244abe50a 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -93,7 +93,7 @@ namespace sys { } // typedefs - typedef std::pair ROI_Key_t; + typedef std::pair ROI_Key_t; typedef std::pair SubROI_Key_t; typedef struct ROIProperties @@ -147,8 +147,8 @@ namespace sys { } TruthProperties_t; // internal containers - std::map< ROI_Key_t,std::vector > ROIMatchedEdepMap; - std::map< ROI_Key_t,std::vector > ROIMatchedHitMap; + std::map< ROI_Key_t, std::vector > ROIMatchedEdepMap; + std::map< ROI_Key_t, std::vector > ROIMatchedHitMap; // some useful functions // geometries @@ -223,7 +223,7 @@ namespace sys { return std::abs(theta); } - // theste are set in the .cc file + // these are set in the .cc file ROIProperties_t CalcROIProperties(recob::Wire const&, size_t const&); std::vector> GetTargetROIs(sim::SimEnergyDeposit const&, double offset); From d9be30652e3ab20f1af8d2c599170400e9440c79 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 1 Apr 2026 16:28:45 -0500 Subject: [PATCH 10/21] Catch and ignore non-active energy depositions. --- sbncode/WireMod/Utility/WireModUtility.cc | 40 +++++++++++++++-------- 1 file changed, 27 insertions(+), 13 deletions(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 3bd0834ed..5fb8e9b94 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -362,6 +362,12 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd { if (edep_ptr->StepLength() == 0) continue; + + geo::TPCID curTPCID; + try { + curTPCID = geometry->PositionToTPC(edep_ptr->MidPoint()).ID(); + } + catch(...) {continue;} // ignore non-active depositions edep_props.x = edep_ptr->X(); edep_props.y = edep_ptr->Y(); @@ -374,8 +380,7 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd total_energy_all += edep_ptr->E(); - const geo::TPCGeo& curTPCGeom = geometry->PositionToTPC(edep_ptr->MidPoint()); - for (auto const& plane : wireReadout->Iterate(curTPCGeom.ID())) { + for (auto const& plane : wireReadout->Iterate(curTPCID)) { int i_p = plane.ID().Plane; auto scales = GetViewScaleValues(edep_props, plane.View()); scales_e_weighted[i_p].r_Q += edep_ptr->E()*scales.r_Q; @@ -456,19 +461,28 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd } edep_col_properties.x_rms_noWeight = std::sqrt(edep_col_properties.x_rms_noWeight); - if (total_energy > 0) + if (total_energy > 0) { edep_col_properties.x_rms = std::sqrt(edep_col_properties.x_rms/total_energy); + } - // projecting edep to plane0, so be mindful when accessing the tick value - const geo::TPCGeo& tpcGeom = geometry->PositionToTPC({edep_col_properties.x, edep_col_properties.y, edep_col_properties.z}); - const auto plane0 = wireReadout->FirstPlane(tpcGeom.ID()); - double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now - edep_col_properties.tick = detPropData.ConvertXToTicks(edep_col_properties.x , plane0.ID()) + offset + tickOffset; - edep_col_properties.tick_rms = ticksPercm*edep_col_properties.x_rms; - edep_col_properties.tick_rms_noWeight = ticksPercm*edep_col_properties.x_rms_noWeight; - edep_col_properties.tick_min = detPropData.ConvertXToTicks(edep_col_properties.x_min, plane0.ID()) + offset + tickOffset; - edep_col_properties.tick_max = detPropData.ConvertXToTicks(edep_col_properties.x_max, plane0.ID()) + offset + tickOffset; - edep_col_properties.total_energy = total_energy; + // get ticks, etc. if the deposition is active + geo::TPCID tpcGeomID; + try { + tpcGeomID = geometry->PositionToTPC({edep_col_properties.x, edep_col_properties.y, edep_col_properties.z}).ID(); + } + catch(...) {} + + if (tpcGeomID.isValid) { + // projecting edep to plane0, so be mindful when accessing the tick value + const auto plane0 = wireReadout->FirstPlane(tpcGeomID); + double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now + edep_col_properties.tick = detPropData.ConvertXToTicks(edep_col_properties.x , plane0.ID()) + offset + tickOffset; + edep_col_properties.tick_rms = ticksPercm*edep_col_properties.x_rms; + edep_col_properties.tick_rms_noWeight = ticksPercm*edep_col_properties.x_rms_noWeight; + edep_col_properties.tick_min = detPropData.ConvertXToTicks(edep_col_properties.x_min, plane0.ID()) + offset + tickOffset; + edep_col_properties.tick_max = detPropData.ConvertXToTicks(edep_col_properties.x_max, plane0.ID()) + offset + tickOffset; + edep_col_properties.total_energy = total_energy; + } return edep_col_properties; From 825d480a48430510406340acf35cc20029adfe67 Mon Sep 17 00:00:00 2001 From: hausnerh Date: Wed, 22 Apr 2026 07:32:54 -0500 Subject: [PATCH 11/21] Fix ThetaZ/angletovert discrepency in calculating the pitch --- sbncode/WireMod/Utility/WireModUtility.hh | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index 244abe50a..e9e147269 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -215,10 +215,11 @@ namespace sys { double ThetaXW(double dxdr, double dydr, double dzdr, double planeAngle) { - double sinPlaneAngle = std::sin(planeAngle); - double cosPlaneAngle = std::cos(planeAngle); + // Want angle to vert, so subtract pi + double sin = std::sin(planeAngle - ::util::pi<>()); + double cos = std::cos(planeAngle - ::util::pi<>()); - double cosG = std::abs(dydr * sinPlaneAngle + dzdr * cosPlaneAngle); + double cosG = std::abs(dydr * sin + dzdr * cos); double theta = std::atan(dxdr / cosG); return std::abs(theta); } From 3e691fb22523f2d53806df6172d7718e4ef47b02 Mon Sep 17 00:00:00 2001 From: hausnerh Date: Tue, 26 May 2026 10:12:06 -0500 Subject: [PATCH 12/21] add pathological check on small q_orig scaling --- sbncode/WireMod/Utility/WireModUtility.cc | 3 +++ 1 file changed, 3 insertions(+) diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 5fb8e9b94..8ca184dc5 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -663,6 +663,9 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, 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; } From 26b359b5959a2cd7594bc41b615ef89353229094 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Tue, 26 May 2026 17:11:42 -0500 Subject: [PATCH 13/21] Hide q_orig warning under verbose. Add in sigma_distance check. --- sbncode/WireMod/Utility/WireModUtility.cc | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 8ca184dc5..fa9f52940 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -630,6 +630,7 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, double q_orig = 0.0; double q_mod = 0.0; double scale_ratio = 1.0; + double sigma_distance = 0.0; // loop over the ticks for(size_t i_t = 0; i_t < roi_data.size(); ++i_t) @@ -638,6 +639,7 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, q_orig = 0.0; q_mod = 0.0; scale_ratio = 1.0; + sigma_distance = 0.0; // loop over the subs for (auto const& subroi_prop : subROIPropVec) @@ -647,12 +649,16 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, 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); + sigma_distance += ((i_t + roi_prop.begin - subroi_prop.center)**2 / subroi_prop.sigma**2)*\ + gausFunc(i_t + roi_prop.begin, subroi_prop.center, subroi_prop.sigma, subroi_prop.total_q); 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; } + sigma_distance = sigma_distance / q_orig; + // do some sanity checks if (isnan(q_orig)) { @@ -664,7 +670,10 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, 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; + if (verbose) std::cout << "WARNING: obtained q_orig < 0.01 ... setting scale to 1" << std::endl; + scale_ratio = 1.0; + } else if (sigma_distance > 3.) { + if (verbose) std::cout << "WARNING: sigma_distance > 3.. setting scale to 1" << std::endl; scale_ratio = 1.0; } else { scale_ratio = q_mod / q_orig; From 9118f6f63e8d6d0f5b0b1f70db79e749e40a77d4 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Tue, 26 May 2026 17:27:09 -0500 Subject: [PATCH 14/21] Fix sq code. --- 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 fa9f52940..bfba735cc 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -649,7 +649,7 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, 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); - sigma_distance += ((i_t + roi_prop.begin - subroi_prop.center)**2 / subroi_prop.sigma**2)*\ + sigma_distance += ((i_t + roi_prop.begin - subroi_prop.center)*(i_t + roi_prop.begin - subroi_prop.center) / subroi_prop.sigma*subroi_prop.sigma)*\ gausFunc(i_t + roi_prop.begin, subroi_prop.center, subroi_prop.sigma, subroi_prop.total_q); if (verbose) From 206a9d0ad9c98e8ffd4c0b50f5b56b834a7147e1 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Tue, 26 May 2026 18:11:48 -0500 Subject: [PATCH 15/21] Fix WM sigma cut to 3^2=9. --- 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 bfba735cc..34b58816d 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -672,7 +672,7 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, } 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; scale_ratio = 1.0; - } else if (sigma_distance > 3.) { + } else if (sigma_distance > 9.) { if (verbose) std::cout << "WARNING: sigma_distance > 3.. setting scale to 1" << std::endl; scale_ratio = 1.0; } else { From 63b9e1cbc755cd073782fd4e51bbc1e0f65b872e Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Tue, 26 May 2026 18:31:03 -0500 Subject: [PATCH 16/21] Fix operator precedence bug. --- 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 34b58816d..138633dba 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -649,7 +649,7 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, 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); - sigma_distance += ((i_t + roi_prop.begin - subroi_prop.center)*(i_t + roi_prop.begin - subroi_prop.center) / subroi_prop.sigma*subroi_prop.sigma)*\ + sigma_distance += ((i_t + roi_prop.begin - subroi_prop.center)*(i_t + roi_prop.begin - subroi_prop.center) / (subroi_prop.sigma*subroi_prop.sigma))*\ gausFunc(i_t + roi_prop.begin, subroi_prop.center, subroi_prop.sigma, subroi_prop.total_q); if (verbose) From 8d37cd1dcb38587b6f251528a2f86ddd20e59073 Mon Sep 17 00:00:00 2001 From: "poppi@baltig.infn.it" Date: Tue, 30 Jun 2026 18:40:33 +0200 Subject: [PATCH 17/21] 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 | 196 ++++++++++++++++--- sbncode/WireMod/Utility/WireModUtility.hh | 28 ++- 4 files changed, 218 insertions(+), 36 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 138633dba..3e58e329c 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) { @@ -614,6 +726,45 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: scales.r_sigma *= splines_Sigma_dEdX [plane]->Eval(truth_props.dedr); } + if (applyXXZAngleScale) + { + if (graph2Ds_Charge_XXZAngle[plane] == nullptr || + 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())); + 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())); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + if (applyXdQdXScale) + { + if (graph2Ds_Charge_XdQdX[plane] == nullptr || + 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); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = graph2Ds_Sigma_XdQdX [plane]->Interpolate(truth_props.x, truth_props.dqdr); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + if (applyXZAngledQdXScale) + { + if (graph2Ds_Charge_XZAngledQdX[plane] == nullptr || + 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); + 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); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + return scales; } @@ -621,11 +772,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; @@ -641,6 +793,8 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, scale_ratio = 1.0; sigma_distance = 0.0; + double t = i_t + roi_prop.begin; + // loop over the subs for (auto const& subroi_prop : subROIPropVec) { @@ -653,45 +807,39 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, gausFunc(i_t + roi_prop.begin, subroi_prop.center, subroi_prop.sigma, subroi_prop.total_q); 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=" << sigma_distance << "sigma" + << (sigma_distance <= sigmaWindow ? " [scaled]" : " [pass-through]") << '\n' + << " q_orig+=" << q_this_orig << '\n'; } - - sigma_distance = sigma_distance / q_orig; + + double delta = q_mod - q_orig; // do some sanity checks - 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 (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 if (verbose) std::cout << "WARNING: obtained q_orig < 0.01 ... setting scale to 1" << std::endl; - scale_ratio = 1.0; } else if (sigma_distance > 9.) { if (verbose) std::cout << "WARNING: sigma_distance > 3.. setting scale to 1" << std::endl; - scale_ratio = 1.0; - } else { - scale_ratio = q_mod / q_orig; + // 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); } - 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; + else { + scale_ratio = q_mod / q_orig; + roi_data[i_t] = scale_ratio * roi_data[i_t]; } - 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; + << ", scale_ratio=" << scale_ratio << " [multiplicative]" + << ", delta=" << delta << " [additive]" << std::endl; } // we're done now diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index e9e147269..ba1d15fa0 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" @@ -43,7 +44,11 @@ namespace sys { bool applyXZAngleScale; // do we scale with XZ angle? bool applyYZAngleScale; // do we scale with YZ angle? bool applydEdXScale; // do we scale with dEdx? + bool applyXXZAngleScale; // do we scale with X vs XZ angle? + bool applyXdQdXScale; // do we scale with X vs dQ/dX? + 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 double readoutWindowTicks; // how many ticks are in the readout window? double tickOffset; // do we want an offset in the ticks? @@ -59,6 +64,13 @@ namespace sys { std::vector splines_Sigma_dEdX; // the splines for the width correction in dEdX std::vector graph2Ds_Charge_YZ; // the graphs for the charge correction in YZ std::vector graph2Ds_Sigma_YZ; // the graphs for the width correction in YZ + + std::vector graph2Ds_Charge_XXZAngle; // the graphs for the charge correction in X vs XZ angle + std::vector graph2Ds_Sigma_XXZAngle; // the graphs for the width correction in X vs XZ angle + std::vector graph2Ds_Charge_XdQdX; // the graphs for charge correction in X vs dQ/dX + std::vector graph2Ds_Sigma_XdQdX; // the graphs for width correction in X vs dQ/dX + std::vector graph2Ds_Charge_XZAngledQdX; // the graphs for charge correction in XZ angle vs dQ/dX + std::vector graph2Ds_Sigma_XZAngledQdX; // the graphs for width correction in XZ angle vs dQ/dX std::vector graph2Ds_Charge_XXW; // the graphs for the charge correction in XXW std::vector graph2Ds_Sigma_XXW; // the graphs for the width correction in XXW @@ -75,6 +87,9 @@ namespace sys { const bool& arg_ApplyXZAngleScale = true, const bool& arg_ApplyYZAngleScale = true, const bool& arg_ApplydEdXScale = true, + const bool& arg_ApplyXXZAngleScale = false, + const bool& arg_ApplyXdQdXScale = false, + const bool& arg_ApplyXZAngledQdXScale = false, const bool& arg_ApplyXXWScale = true, const double& arg_TickOffset = 0) : geometry(geom), @@ -86,7 +101,11 @@ namespace sys { applyXZAngleScale(arg_ApplyXZAngleScale), applyYZAngleScale(arg_ApplyYZAngleScale), applydEdXScale(arg_ApplydEdXScale), + applyXXZAngleScale(arg_ApplyXXZAngleScale), + applyXdQdXScale(arg_ApplyXdQdXScale), + applyXZAngledQdXScale(arg_ApplyXZAngledQdXScale), 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 { @@ -142,6 +161,7 @@ namespace sys { double dxdr; double dydr; double dzdr; + double dqdr; double dedr; ScaleValues_t scales_avg[3]; } TruthProperties_t; @@ -200,6 +220,7 @@ namespace sys { double theta = std::atan2(dydrPlaneRel, dzdrPlaneRel); return FoldAngle(theta); } + // theste are set in the .cc file double ThetaXY_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) { @@ -226,6 +247,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 +256,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 +272,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 645ea905f0300383a8a3b5253d1f5cdf48ead960 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Tue, 30 Jun 2026 17:03:06 -0500 Subject: [PATCH 18/21] 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 bb09f01200f1bfaf2587ee1f4965a0c58b3a9d32 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 1 Jul 2026 08:57:04 -0500 Subject: [PATCH 19/21] 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 8c006e8163a69833aa3a5c7c3ee8cef7bdc49324 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Thu, 2 Jul 2026 14:08:44 -0500 Subject: [PATCH 20/21] 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 78c362b149e9a5dfc5b14c3158838926928b760e Mon Sep 17 00:00:00 2001 From: "poppi@baltig.infn.it" Date: Wed, 8 Jul 2026 23:53:59 +0200 Subject: [PATCH 21/21] 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 {