From 382285a626b34886c763492fed2904bdc43eeb56 Mon Sep 17 00:00:00 2001 From: Harry Hausner Date: Wed, 5 Jun 2024 14:10:18 -0500 Subject: [PATCH 1/6] 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 9b7cfc5ed..1c0208f79 100644 --- a/sbncode/CMakeLists.txt +++ b/sbncode/CMakeLists.txt @@ -20,6 +20,7 @@ add_subdirectory(GeometryTools) add_subdirectory(CosmicID) add_subdirectory(DetSim) add_subdirectory(Cluster3D) +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 d8db1d8b0d0f3787ef17184120be8df478351837 Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Wed, 19 Mar 2025 13:59:33 -0500 Subject: [PATCH 2/6] 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 f37272b9f4c193d1bb80329876d48a88cff2b3f8 Mon Sep 17 00:00:00 2001 From: kjplows Date: Wed, 30 Jul 2025 15:36:57 -0500 Subject: [PATCH 3/6] Update dependencies --- CMakeLists.txt | 2 +- ups/product_deps | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 176b4d3c9..6d2f8b3ef 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,7 +16,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR) find_package(cetmodules 3.20.00 REQUIRED) -project(sbncode VERSION 10.06.00.03 LANGUAGES CXX) +project(sbncode VERSION 10.06.00.04 LANGUAGES CXX) message(STATUS "\n\n ========================== ${PROJECT_NAME} ==========================") diff --git a/ups/product_deps b/ups/product_deps index efa661c7b..4d5dd8cb9 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -254,7 +254,7 @@ libdir fq_dir lib product version qual flags genie_xsec v3_04_00 - larcv2 v2_2_6 - -larsoft v10_06_00_01 - +larsoft v10_06_00_02 - sbnalg v10_06_00_03 - sbndaq_artdaq_core v1_10_06 - sbndata v01_07 - From 9e8d0366543c5759b98b45abd24cbd9672f095fc Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Thu, 16 Oct 2025 17:36:38 -0500 Subject: [PATCH 4/6] WireMod for sbnd production branch --- sbncode/CMakeLists.txt | 2 + sbncode/WireMod/CMakeLists.txt | 1 + sbncode/WireMod/Utility/CMakeLists.txt | 45 ++ sbncode/WireMod/Utility/WireModUtility.cc | 654 ++++++++++++++++++++++ sbncode/WireMod/Utility/WireModUtility.hh | 228 ++++++++ 5 files changed, 930 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..6b2c41557 100644 --- a/sbncode/CMakeLists.txt +++ b/sbncode/CMakeLists.txt @@ -48,3 +48,5 @@ execute_process(COMMAND sh ${CMAKE_CURRENT_SOURCE_DIR}/Supera/setup.sh icarus add_subdirectory(Supera) add_subdirectory(LArG4) + +add_subdirectory(WireMod) 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..89b4440f4 --- /dev/null +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -0,0 +1,654 @@ +#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 (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 + 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; +} + +//--- 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) +{ + // 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()); + 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; + + // 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 (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; + } + } + + 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}); + 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; + + + 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 + auto const& plane_obj = wireReadout->Plane(curTPCGeomPtr->ID(), view); + unsigned int plane = plane_obj.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, 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) + { + 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, 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) + { + 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..9fabdb1b2 --- /dev/null +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -0,0 +1,228 @@ +#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 "larcorealg/Geometry/WireReadoutGeom.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 "lardataobj/Simulation/SimChannel.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 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? + 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 geo::WireReadoutGeom* wireRead, + 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), + wireReadout(wireRead), + 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 + // 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, const geo::PlaneGeo& plane, const geo::TPCGeo& tpcGeom, double offset = 0) + { + double tick = planeXToTick(xPos, plane, 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&); + void FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, double offset); + + 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 75853f7ae017cda77b80d0be0f7a65b48b43e101 Mon Sep 17 00:00:00 2001 From: hausnerh Date: Wed, 22 Oct 2025 14:36:16 -0500 Subject: [PATCH 5/6] Update WireMod with more possible ways to scale --- sbncode/WireMod/Utility/WireModUtility.cc | 49 ++++++++++++++++++++++- sbncode/WireMod/Utility/WireModUtility.hh | 17 ++++++++ 2 files changed, 65 insertions(+), 1 deletion(-) diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index d61bac5ac..35e54ecc9 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -361,6 +361,8 @@ 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.dqdr = edep_ptr->NumElectrons() / edep_ptr->StepLength(); edep_props.dedr = edep_ptr->E() / edep_ptr->StepLength(); total_energy_all += edep_ptr->E(); @@ -408,6 +410,8 @@ 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.dqdr = 0.; edep_col_properties.dedr = 0.; double total_energy = 0.0; @@ -426,6 +430,8 @@ 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.dqdr += edep_ptr->E()*edep_ptr->NumElectrons() / edep_ptr->StepLength(); edep_col_properties.dedr += edep_ptr->E()*edep_ptr->E() / edep_ptr->StepLength(); } @@ -437,6 +443,8 @@ 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.dqdr = edep_col_properties.dqdr / total_energy; edep_col_properties.dedr = edep_col_properties.dedr / total_energy; } @@ -566,7 +574,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 +584,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..1ae74c152 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -42,6 +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? + 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? double readoutWindowTicks; // how many ticks are in the readout window? double tickOffset; // do we want an offset in the ticks? @@ -58,6 +61,13 @@ 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 + 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 + // 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 +81,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 double& arg_TickOffset = 0) : geometry(geom), wireReadout(wireRead), @@ -81,6 +94,9 @@ namespace sys { applyXZAngleScale(arg_ApplyXZAngleScale), applyYZAngleScale(arg_ApplyYZAngleScale), applydEdXScale(arg_ApplydEdXScale), + applyXXZAngleScale(arg_ApplyXXZAngleScale), + applyXdQdXScale(arg_ApplyXdQdXScale), + applyXZAngledQdXScale(arg_ApplyXZAngledQdXScale), 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 +152,7 @@ namespace sys { double dxdr; double dydr; double dzdr; + double dqdr; double dedr; ScaleValues_t scales_avg[3]; } TruthProperties_t; From d3fd5db25021d6de3715257a48819a7a482874c9 Mon Sep 17 00:00:00 2001 From: Thomas Wester Date: Thu, 23 Oct 2025 12:01:41 -0500 Subject: [PATCH 6/6] fix build --- sbncode/CMakeLists.txt | 2 -- sbncode/WireMod/Utility/WireModUtility.hh | 15 --------------- 2 files changed, 17 deletions(-) diff --git a/sbncode/CMakeLists.txt b/sbncode/CMakeLists.txt index 3440f348d..14a8e0265 100644 --- a/sbncode/CMakeLists.txt +++ b/sbncode/CMakeLists.txt @@ -49,5 +49,3 @@ execute_process(COMMAND sh ${CMAKE_CURRENT_SOURCE_DIR}/Supera/setup.sh icarus add_subdirectory(Supera) add_subdirectory(LArG4) - -add_subdirectory(WireMod) 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;