From 382285a626b34886c763492fed2904bdc43eeb56 Mon Sep 17 00:00:00 2001 From: Harry Hausner Date: Wed, 5 Jun 2024 14:10:18 -0500 Subject: [PATCH] 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