diff --git a/sbncode/CAFMaker/RecoUtils/CMakeLists.txt b/sbncode/CAFMaker/RecoUtils/CMakeLists.txt index b6abdbdb9..5401c9bbd 100644 --- a/sbncode/CAFMaker/RecoUtils/CMakeLists.txt +++ b/sbncode/CAFMaker/RecoUtils/CMakeLists.txt @@ -2,6 +2,7 @@ art_make_library( LIBRARY_NAME caf_RecoUtils SOURCE RecoUtils.cc LIBRARIES + sbnanaobj::StandardRecord art::Framework_Core art::Framework_Services_Registry art::Framework_Principal diff --git a/sbncode/CMakeLists.txt b/sbncode/CMakeLists.txt index daed2590a..8b9a8ff9f 100644 --- a/sbncode/CMakeLists.txt +++ b/sbncode/CMakeLists.txt @@ -21,6 +21,7 @@ add_subdirectory(CosmicID) add_subdirectory(DetSim) add_subdirectory(Cluster3D) add_subdirectory(HitFinder) +add_subdirectory(WireMod) # Supera # diff --git a/sbncode/HitFinder/ChannelROIToWire_module.cc b/sbncode/HitFinder/ChannelROIToWire_module.cc index 42d304537..497ffedcb 100644 --- a/sbncode/HitFinder/ChannelROIToWire_module.cc +++ b/sbncode/HitFinder/ChannelROIToWire_module.cc @@ -51,6 +51,7 @@ class ChannelROIToWire : public art::EDProducer std::vector fWireModuleLabelVec; ///< vector of modules that made digits std::vector fOutInstanceLabelVec; ///< The output instance labels to apply bool fDiagnosticOutput; ///< secret diagnostics flag + bool fUseFloatPrecision; ///< use SignalROIF() to preserve sub-ADC precision (no rounding) size_t fEventCount; ///< count of event processed const geo::WireReadoutGeom* fChannelMapAlg = &art::ServiceHandle()->Get(); @@ -77,6 +78,7 @@ void ChannelROIToWire::reconfigure(fhicl::ParameterSet const& pset) fWireModuleLabelVec = pset.get>("WireModuleLabelVec", std::vector()={"decon1droi"}); fOutInstanceLabelVec = pset.get> ("OutInstanceLabelVec", {"PHYSCRATEDATA"}); fDiagnosticOutput = pset.get< bool >("DiagnosticOutput", false); + fUseFloatPrecision = pset.get< bool >("UseFloatPrecision", false); if (fWireModuleLabelVec.size() != fOutInstanceLabelVec.size()) { @@ -131,19 +133,24 @@ void ChannelROIToWire::produce(art::Event& evt) // Create an ROI vector for output recob::Wire::RegionsOfInterest_t ROIVec; - - // Loop through the ROIs for this channel - const recob::ChannelROI::RegionsOfInterest_t& channelROIs = channelROI.SignalROI(); - for(const auto& range : channelROIs.get_ranges()) + if (fUseFloatPrecision) { - size_t startTick = range.begin_index(); - - std::vector dataVec(range.data().size()); - - for(size_t binIdx = 0; binIdx < range.data().size(); binIdx++) dataVec[binIdx] = std::round(range.data()[binIdx] / ADCScaleFactor); - - ROIVec.add_range(startTick, std::move(dataVec)); + const recob::ChannelROI::RegionsOfInterest_f& channelROIsF = channelROI.SignalROIF(); + for(const auto& range : channelROIsF.get_ranges()) + ROIVec.add_range(range.begin_index(), range.data()); + } + else + { + const recob::ChannelROI::RegionsOfInterest_t& channelROIs = channelROI.SignalROI(); + for(const auto& range : channelROIs.get_ranges()) + { + size_t startTick = range.begin_index(); + std::vector dataVec(range.data().size()); + for(size_t binIdx = 0; binIdx < range.data().size(); binIdx++) + dataVec[binIdx] = std::round(range.data()[binIdx] / ADCScaleFactor); + ROIVec.add_range(startTick, std::move(dataVec)); + } } wireCol->push_back(recob::WireCreator(std::move(ROIVec),channel,view).move()); diff --git a/sbncode/WireMod/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..42a65ae0d --- /dev/null +++ b/sbncode/WireMod/Utility/CMakeLists.txt @@ -0,0 +1,49 @@ +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 + lardataobj::Simulation + nusimdata::SimulationBase + larevt::SpaceCharge + sbnobj::ICARUS_TPC + 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..001ef7ba0 --- /dev/null +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -0,0 +1,1168 @@ +#include "sbncode/WireMod/Utility/WireModUtility.hh" +#include "lardataalg/DetectorInfo/DetectorPropertiesData.h" + +#include +#include +#include +#include "TVector3.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_tfloat manually. + auto const& roi_short = chanROI.SignalROI().get_ranges()[roi_idx]; + const float scale = 1.0f / static_cast(chanROI.ADCScaleFactor()); + + ROIProperties_t roi_vals; + roi_vals.channel = chanROI.Channel(); + roi_vals.view = wireReadout->View(chanROI.Channel()); + roi_vals.begin = roi_short.begin_index(); + roi_vals.end = roi_short.end_index(); + roi_vals.center = 0; + roi_vals.total_q = 0; + roi_vals.sigma = 0; + + auto const& short_data = roi_short.data(); + for (size_t i_t = 0; i_t < short_data.size(); ++i_t) + { + float val = static_cast(short_data[i_t]) * scale; + roi_vals.center += val * (i_t + roi_vals.begin); + roi_vals.total_q += val; + } + roi_vals.center = roi_vals.center / roi_vals.total_q; + + for (size_t i_t = 0; i_t < short_data.size(); ++i_t) + { + float val = static_cast(short_data[i_t]) * scale; + roi_vals.sigma += val * (i_t + roi_vals.begin - roi_vals.center) + * (i_t + roi_vals.begin - roi_vals.center); + } + roi_vals.sigma = std::sqrt(roi_vals.sigma / roi_vals.total_q); + if (roi_vals.end - roi_vals.begin == 1) + { + roi_vals.center += 0.5; + roi_vals.sigma = 0.5; + } + + return roi_vals; +} + +//--- GetTargetROIs --- +std::vector> sys::WireModUtility::GetTargetROIs(sim::SimEnergyDeposit const& shifted_edep, double offset) +{ + // 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 --- +// SimChannel/IDE analogue of FillROIMatchedEdepMap. Unlike an edep, a sim::IDE already +// carries its readout channel (SimChannel::Channel()) and time (the TDCIDEMap TDC), so no +// geometric X->tick projection is needed -- the TDC is converted to a readout tick and the +// covering ROI is found directly. Builds the flat fIDEVec and the ROI->index map. +void sys::WireModUtility::FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, detinfo::DetectorClocksData const& clockData, double offset) +{ + ROIMatchedIDEMap.clear(); + fIDEVec.clear(); + + // map channel -> index into wireVec + std::unordered_map wireChannelMap; + for (size_t i_w = 0; i_w < wireVec.size(); ++i_w) + wireChannelMap[wireVec[i_w].Channel()] = i_w; + + for (auto const& simch : simchVec) + { + raw::ChannelID_t channel = simch.Channel(); + auto wcIt = wireChannelMap.find(channel); + if (wcIt == wireChannelMap.end()) + continue; + auto const& target_wire = wireVec.at(wcIt->second); + + // wire closest to the charge (used later for the SCE-undo) + geo::WireID wireID; + auto const wires = wireReadout->ChannelToWire(channel); + if (!wires.empty()) wireID = wires[0]; + + for (auto const& tdcide : simch.TDCIDEMap()) + { + // convert the SimChannel TDC to a readout tick (same frame as the ROI samples) + double tick = clockData.TPCTDC2Tick(tdcide.first) + offset + tickOffset; + if (tick < 0) continue; + unsigned int sample = (unsigned int) std::round(tick); + + if (not target_wire.SignalROI().is_valid() || + target_wire.SignalROI().empty() || + target_wire.SignalROI().n_ranges() == 0 || + target_wire.SignalROI().size() <= sample || + target_wire.SignalROI().is_void(sample) ) + continue; + + auto range_number = target_wire.SignalROI().find_range_iterator(sample) - target_wire.SignalROI().begin_range(); + ROI_Key_t roi_key = std::make_pair(target_wire.Channel(), range_number); + + for (sim::IDE const& ide : tdcide.second) + { + fIDEVec.push_back( MatchedIDE_t{channel, wireID, tick, &ide} ); + ROIMatchedIDEMap[roi_key].push_back(fIDEVec.size() - 1); + } + } + } +} + +//--- FillROIMatchedIDEMap (ChannelROI overload) --- +void sys::WireModUtility::FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& chanROIVec, detinfo::DetectorClocksData const& clockData, double offset) +{ + ROIMatchedIDEMap.clear(); + fIDEVec.clear(); + + std::unordered_map chanROIChannelMap; + for (size_t i_c = 0; i_c < chanROIVec.size(); ++i_c) + chanROIChannelMap[chanROIVec[i_c].Channel()] = i_c; + + for (auto const& simch : simchVec) + { + raw::ChannelID_t channel = simch.Channel(); + auto ccIt = chanROIChannelMap.find(channel); + if (ccIt == chanROIChannelMap.end()) + continue; + auto const& target_chanROI = chanROIVec.at(ccIt->second); + auto const& roiShort = target_chanROI.SignalROI(); + + geo::WireID wireID; + auto const wires = wireReadout->ChannelToWire(channel); + if (!wires.empty()) wireID = wires[0]; + + for (auto const& tdcide : simch.TDCIDEMap()) + { + double tick = clockData.TPCTDC2Tick(tdcide.first) + offset + tickOffset; + if (tick < 0) continue; + unsigned int sample = (unsigned int) std::round(tick); + + if (not roiShort.is_valid() || + roiShort.empty() || + roiShort.n_ranges() == 0 || + roiShort.size() <= sample || + roiShort.is_void(sample) ) + continue; + + auto range_number = roiShort.find_range_iterator(sample) - roiShort.begin_range(); + ROI_Key_t roi_key = std::make_pair(target_chanROI.Channel(), range_number); + + for (sim::IDE const& ide : tdcide.second) + { + fIDEVec.push_back( MatchedIDE_t{channel, wireID, tick, &ide} ); + ROIMatchedIDEMap[roi_key].push_back(fIDEVec.size() - 1); + } + } + } +} + +//--- FillROIMatchedEdepMap --- +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); + } + } +} + +//--- FillROIMatchedEdepMap (ChannelROI overload) --- +void sys::WireModUtility::FillROIMatchedEdepMap(std::vector const& edepVec, std::vector const& chanROIVec, double offset) +{ + ROIMatchedEdepMap.clear(); + + std::unordered_map chanROIChannelMap; + for (size_t i_c = 0; i_c < chanROIVec.size(); ++i_c) + chanROIChannelMap[chanROIVec[i_c].Channel()] = i_c; + + for (size_t i_e = 0; i_e < edepVec.size(); ++i_e) + { + std::vector> target_rois = GetTargetROIs(edepVec[i_e], offset); + + for (auto const& target_roi : target_rois) + { + if (chanROIChannelMap.find(target_roi.first) == chanROIChannelMap.end()) + continue; + + auto const& target_chanROI = chanROIVec.at(chanROIChannelMap[target_roi.first]); + + // SignalROIF() returns by value, so use SignalROI() instead. + auto const& roiShort = target_chanROI.SignalROI(); + if (not roiShort.is_valid() || + roiShort.empty() || + roiShort.n_ranges() == 0 || + roiShort.size() <= target_roi.second || + roiShort.is_void(target_roi.second) ) + continue; + + auto range_number = roiShort.find_range_iterator(target_roi.second) - roiShort.begin_range(); + + ROIMatchedEdepMap[std::make_pair(target_chanROI.Channel(), range_number)].push_back(i_e); + } + } +} + +//--- FillROIMatchedHitMap (ChannelROI overload) --- +void sys::WireModUtility::FillROIMatchedHitMap(std::vector const& hitVec, std::vector const& chanROIVec) +{ + ROIMatchedHitMap.clear(); + + std::unordered_map chanROIChannelMap; + for (size_t i_c = 0; i_c < chanROIVec.size(); ++i_c) + chanROIChannelMap[chanROIVec[i_c].Channel()] = i_c; + + for (size_t i_h = 0; i_h < hitVec.size(); ++i_h) + { + std::vector> target_rois = GetHitTargetROIs(hitVec[i_h]); + + for (auto const& target_roi : target_rois) + { + if (chanROIChannelMap.find(target_roi.first) == chanROIChannelMap.end()) + continue; + + auto const& target_chanROI = chanROIVec.at(chanROIChannelMap[target_roi.first]); + + // SignalROIF() returns by value, so use SignalROI() instead. + auto const& roiShort = target_chanROI.SignalROI(); + if (not roiShort.is_valid() || + roiShort.empty() || + roiShort.n_ranges() == 0 || + roiShort.size() <= target_roi.second || + roiShort.is_void(target_roi.second) ) + continue; + + auto range_number = roiShort.find_range_iterator(target_roi.second) - roiShort.begin_range(); + + ROIMatchedHitMap[std::make_pair(target_chanROI.Channel(), range_number)].push_back(i_h); + } + } +} + +//--- CalcSubROIProperties --- +std::vector sys::WireModUtility::CalcSubROIProperties(sys::WireModUtility::ROIProperties_t const& roi_properties, std::vector const& hitPtrVec) +{ + 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(); + } + + // 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()); + std::map planeProjEdepTick; + for (const auto& planeGeom : wireReadout->Iterate(curTPCGeom.ID())) + { + double projTick = detPropData.ConvertXToTicks(edep_ptr->X(), planeGeom.ID()) + offset + tickOffset; + readout::ROPID ropID = wireReadout->WirePlaneToROP(planeGeom.ID()); + if (planeProjEdepTick.count(ropID) == 1) + { + std::cout << "MatchEdepsToSubROIs - Warning, overwriting tick projection at" << ropID.toString() + << "\n Changing from " << planeProjEdepTick[ropID] << " to " << projTick << std::endl; + } + planeProjEdepTick[ropID] = projTick; + } + + // loop over subROIs + unsigned int closest_hit = std::numeric_limits::max(); + 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]; + readout::ROPID subroiROP = wireReadout->ChannelToROP(subroi_prop.channel); + double edep_tick = planeProjEdepTick.at(subroiROP); + if (edep_tick > subroi_prop.center-subroi_prop.sigma && edep_tick < subroi_prop.center+subroi_prop.sigma) + { + EDepMatchedSubROIMap[i_e].push_back(i_h); + 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; + + geo::TPCID curTPCID; + try { + curTPCID = geometry->PositionToTPC(edep_ptr->MidPoint()).ID(); + } + catch(...) {continue;} // ignore non-active depositions + + edep_props.x = edep_ptr->X(); + edep_props.y = edep_ptr->Y(); + 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(); + + for (auto const& plane : wireReadout->Iterate(curTPCID)) { + int i_p = plane.ID().Plane; + auto scales = GetViewScaleValues(edep_props, plane.View()); + scales_e_weighted[i_p].r_Q += edep_ptr->E()*scales.r_Q; + 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); + } + + // get ticks, etc. if the deposition is active + geo::TPCID tpcGeomID; + try { + tpcGeomID = geometry->PositionToTPC({edep_col_properties.x, edep_col_properties.y, edep_col_properties.z}).ID(); + } + catch(...) {} + + if (tpcGeomID.isValid) { + // projecting edep to plane0, so be mindful when accessing the tick value + const auto plane0 = wireReadout->FirstPlane(tpcGeomID); + double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now + edep_col_properties.tick = detPropData.ConvertXToTicks(edep_col_properties.x , plane0.ID()) + offset + tickOffset; + edep_col_properties.tick_rms = ticksPercm*edep_col_properties.x_rms; + edep_col_properties.tick_rms_noWeight = ticksPercm*edep_col_properties.x_rms_noWeight; + edep_col_properties.tick_min = detPropData.ConvertXToTicks(edep_col_properties.x_min, plane0.ID()) + offset + tickOffset; + edep_col_properties.tick_max = detPropData.ConvertXToTicks(edep_col_properties.x_max, plane0.ID()) + offset + tickOffset; + edep_col_properties.total_energy = total_energy; + } + + + return edep_col_properties; +} + +//--- WireToTrajectoryPosition (SCE: at-the-wire -> true trajectory frame) --- +geo::Point_t sys::WireModUtility::WireToTrajectoryPosition(const geo::Point_t& loc, const geo::TPCID& tpc) const +{ + geo::Point_t ret = loc; + if (fSCE && fSCE->EnableSimSpatialSCE()) + { + geo::Vector_t offset = fSCE->GetCalPosOffsets(ret, tpc.TPC); + ret.SetX(ret.X() + offset.X()); + ret.SetY(ret.Y() + offset.Y()); + ret.SetZ(ret.Z() + offset.Z()); + } + return ret; +} + +//--- TrajectoryToWirePosition (SCE: true -> at-the-wire frame) --- +geo::Point_t sys::WireModUtility::TrajectoryToWirePosition(const geo::Point_t& loc, const geo::Vector_t& driftdir) const +{ + geo::Point_t ret = loc; + int corr = driftdir.X(); // returned X offset is the drift; sign it by the drift direction + if (fSCE && fSCE->EnableSimSpatialSCE()) + { + geo::Vector_t offset = fSCE->GetPosOffsets(ret); + ret.SetX(ret.X() + corr * offset.X()); + ret.SetY(ret.Y() + offset.Y()); + ret.SetZ(ret.Z() + offset.Z()); + } + return ret; +} + +//--- MatchIDEsToSubROIs --- +// SimChannel/IDE analogue of MatchEdepsToSubROIs. Each IDE already carries its readout tick +// (no per-plane X->tick projection needed), so it is assigned to the closest sub-ROI whose +// distance |tick - center|/sigma is below the same threshold used for edeps (5 sigma). +std::map> +sys::WireModUtility::MatchIDEsToSubROIs(std::vector const& subROIPropVec, + std::vector const& idePtrVec) +{ + std::map> SubROIMatchedIDEMap; // subROI idx -> ide idxs + + for (unsigned int i_e = 0; i_e < idePtrVec.size(); ++i_e) + { + double ide_tick = idePtrVec[i_e]->tick; + unsigned int closest = std::numeric_limits::max(); + float min_dist = std::numeric_limits::max(); + for (unsigned int i_h = 0; i_h < subROIPropVec.size(); ++i_h) + { + auto const& s = subROIPropVec[i_h]; + if (s.sigma <= 0) continue; + float dist = std::abs((float)(ide_tick - s.center)) / s.sigma; + if (dist < min_dist) { min_dist = dist; closest = i_h; } + } + if (closest != std::numeric_limits::max() && min_dist < 5) + SubROIMatchedIDEMap[closest].push_back(i_e); + } + + std::map> ReturnMap; + for (auto const& kv : SubROIMatchedIDEMap) + { + auto key = subROIPropVec[kv.first].key; + for (auto const& i_e : kv.second) + ReturnMap[key].push_back(idePtrVec[i_e]); + } + return ReturnMap; +} + +//--- CalcPropertiesFromIDEs --- +// SimChannel/IDE analogue of CalcPropertiesFromEdeps. Picks the dominant track (max summed +// energy, abs(trackID) as in TrackCaloSkimmer), builds the numElectrons-weighted centroid in +// the SCE-undone (trajectory) frame, and reads the direction (dxdr/dydr/dzdr) from the matched +// MCParticle's nearest trajectory point. dedr/dqdr use the trajectory-derived pitch. The +// scales_avg pre-pass of the edep version is dropped (those fields are never read downstream). +sys::WireModUtility::TruthProperties_t +sys::WireModUtility::CalcPropertiesFromIDEs(std::vector const& idePtrVec, + std::map const& particleMap) +{ + TruthProperties_t props; + props.x = props.y = props.z = 0.f; + props.x_rms = props.x_rms_noWeight = 0.f; + props.x_min = std::numeric_limits::max(); + props.x_max = std::numeric_limits::lowest(); + props.tick = props.tick_rms = props.tick_rms_noWeight = 0.f; + props.tick_min = std::numeric_limits::max(); + props.tick_max = std::numeric_limits::lowest(); + props.dxdr = props.dydr = props.dzdr = 0.; + props.dedr = props.dqdr = 0.; + props.total_energy = 0.f; + for (size_t i_p = 0; i_p < 3; ++i_p) { props.scales_avg[i_p].r_Q = 1.; props.scales_avg[i_p].r_sigma = 1.; } + if (idePtrVec.empty()) return props; + + // dominant track = max summed energy (abs(trackID), matching TrackCaloSkimmer) + std::map energy_per_trkid; + for (auto const& m : idePtrVec) + energy_per_trkid[std::abs(m->ide->trackID)] += m->ide->energy; + int trkid_max = 0; double energy_max = -1.; + for (auto const& e : energy_per_trkid) + if (e.second > energy_max) { energy_max = e.second; trkid_max = e.first; } + + // numElectrons-weighted centroid in the SCE-undone frame + energy/electron/tick sums + struct Acc { double xtraj, tick, ne; }; + std::vector accs; + double sumNe = 0., sumE = 0., cx = 0., cy = 0., cz = 0., tsum = 0.; + geo::WireID domWire; + for (auto const& m : idePtrVec) + { + if (std::abs(m->ide->trackID) != trkid_max) continue; + double ne = m->ide->numElectrons; + double e = m->ide->energy; + geo::Point_t p_raw(m->ide->x, m->ide->y, m->ide->z); + geo::Point_t p_traj = WireToTrajectoryPosition(p_raw, m->wire); + cx += ne * p_traj.X(); cy += ne * p_traj.Y(); cz += ne * p_traj.Z(); + tsum += ne * m->tick; + sumNe += ne; sumE += e; + accs.push_back({p_traj.X(), m->tick, ne}); + if (p_traj.X() < props.x_min) props.x_min = p_traj.X(); + if (p_traj.X() > props.x_max) props.x_max = p_traj.X(); + if (m->tick < props.tick_min) props.tick_min = m->tick; + if (m->tick > props.tick_max) props.tick_max = m->tick; + domWire = m->wire; + } + if (sumNe <= 0) return props; + cx /= sumNe; cy /= sumNe; cz /= sumNe; + props.x = cx; props.y = cy; props.z = cz; + props.total_energy = sumE; + props.tick = tsum / sumNe; + + double xr = 0., xrnw = 0., tr = 0.; + for (auto const& a : accs) + { + xr += a.ne * (a.xtraj - cx) * (a.xtraj - cx); + xrnw += (a.xtraj - cx) * (a.xtraj - cx); + tr += a.ne * (a.tick - props.tick) * (a.tick - props.tick); + } + props.x_rms = std::sqrt(xr / sumNe); + props.x_rms_noWeight = std::sqrt(xrnw); + props.tick_rms = std::sqrt(tr / sumNe); + props.tick_rms_noWeight = props.x_rms_noWeight; + + // direction / pitch from the matched MCParticle trajectory (nearest point to centroid) + auto pit = particleMap.find(trkid_max); + if (pit != particleMap.end() && pit->second != nullptr) + { + const simb::MCParticle* part = pit->second; + TVector3 hp(cx, cy, cz); + double closest = -1.; TVector3 dir; + for (unsigned int i = 0; i < part->NumberTrajectoryPoints(); ++i) + { + double d = (part->Position(i).Vect() - hp).Mag(); + if (closest < 0. || d < closest) + { + closest = d; + dir = part->Momentum(i).Vect().Unit(); + } + } + if (closest >= 0. && dir.Mag() > 1e-4) + { + props.dxdr = dir.X(); + props.dydr = dir.Y(); + props.dzdr = dir.Z(); + + if (domWire.isValid) + { + geo::PlaneID plane(domWire.Cryostat, domWire.TPC, domWire.Plane); + geo::PlaneGeo const& pg = wireReadout->Plane(plane); + double angletovert = wireReadout->WireAngleToVertical(pg.View(), plane) - 0.5 * ::util::pi<>(); + double cosgamma = std::abs(std::cos(angletovert) * dir.Z() + std::sin(angletovert) * dir.Y()); + double pitch = (cosgamma > 1e-6) ? pg.WirePitch() / cosgamma : pg.WirePitch(); + if (pitch > 0) + { + props.dedr = sumE / pitch; // dE/dx [MeV/cm] + props.dqdr = sumNe / pitch; // dQ/dx [electrons/cm] + } + } + } + } + + return props; +} + +// TGraph2D lookup: interpolate=true uses Delaunay (clamped), false uses nearest bin center. +static double g2dLookup(TGraph2D* g, double x, double y, bool interpolate) +{ + if (interpolate) + return g->Interpolate(std::clamp(x, g->GetXmin(), g->GetXmax()), + std::clamp(y, g->GetYmin(), g->GetYmax())); + + int n = g->GetN(); + double* xs = g->GetX(); + double* ys = g->GetY(); + double* zs = g->GetZ(); + + std::vector ux(xs, xs + n), uy(ys, ys + n); + std::sort(ux.begin(), ux.end()); + ux.erase(std::unique(ux.begin(), ux.end()), ux.end()); + std::sort(uy.begin(), uy.end()); + uy.erase(std::unique(uy.begin(), uy.end()), uy.end()); + + double qx = std::clamp(x, ux.front(), ux.back()); + double qy = std::clamp(y, uy.front(), uy.back()); + + auto findBin = [](std::vector const& c, double q) -> int { + int m = static_cast(c.size()); + if (m == 1) return 0; + double left = c[0] - (c[1] - c[0]) * 0.5; + for (int i = 0; i < m; ++i) { + double right = 2.0 * c[i] - left; + if (q < right || i == m - 1) return i; + left = right; + } + return m - 1; + }; + + double cx = ux[findBin(ux, qx)]; + double cy = uy[findBin(uy, qy)]; + + for (int i = 0; i < n; ++i) + if (std::abs(xs[i] - cx) < 1e-9 && std::abs(ys[i] - cy) < 1e-9) + return zs[i]; + + return 1.0; // bin not populated in the map (sparse grid) — neutral scale +} + +//--- GetScaleValues --- +sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetScaleValues(sys::WireModUtility::TruthProperties_t const& truth_props, sys::WireModUtility::ROIProperties_t const& roi_vals) +{ + 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 = g2dLookup(graph2Ds_Charge_YZ[plane], truth_props.z, truth_props.y, useGraph2DInterpolation); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = g2dLookup(graph2Ds_Sigma_YZ [plane], truth_props.z, truth_props.y, useGraph2DInterpolation); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + if (applyXXWScale) + { + if (graph2Ds_Charge_XXW[plane] == nullptr || + graph2Ds_Sigma_XXW [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply XXW scale factor, but could not find graphs. Check that you have set those in the utility."; + double thXW = ThetaXW(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()); + temp_scale = g2dLookup(graph2Ds_Charge_XXW[plane], truth_props.x, thXW, useGraph2DInterpolation); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = g2dLookup(graph2Ds_Sigma_XXW [plane], truth_props.x, thXW, useGraph2DInterpolation); + 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); + } + + 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 = g2dLookup(graph2Ds_Charge_XXZAngle[plane], truth_props.x, ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), useGraph2DInterpolation); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = g2dLookup(graph2Ds_Sigma_XXZAngle [plane], truth_props.x, ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), useGraph2DInterpolation); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + 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 = g2dLookup(graph2Ds_Charge_XdQdX[plane], truth_props.x, truth_props.dqdr, useGraph2DInterpolation); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = g2dLookup(graph2Ds_Sigma_XdQdX [plane], truth_props.x, truth_props.dqdr, useGraph2DInterpolation); + 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 = g2dLookup(graph2Ds_Charge_XZAngledQdX[plane], ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), truth_props.dqdr, useGraph2DInterpolation); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = g2dLookup(graph2Ds_Sigma_XZAngledQdX [plane], ThetaXZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()), truth_props.dqdr, useGraph2DInterpolation); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + + 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, + double sigmaWindow) +{ + // do you want a bunch of messages? + bool verbose = false; + + // initialize some values + double q_orig = 0.0; + double q_mod = 0.0; + double scale_ratio = 1.0; + double sigma_distance = 0.0; + + // loop over the ticks + for(size_t i_t = 0; i_t < roi_data.size(); ++i_t) + { + // reset your values + q_orig = 0.0; + q_mod = 0.0; + scale_ratio = 1.0; + sigma_distance = 0.0; + + double t = i_t + roi_prop.begin; + + // loop over the subs + for (auto const& subroi_prop : subROIPropVec) + { + // 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); + sigma_distance += ((i_t + roi_prop.begin - subroi_prop.center)*(i_t + roi_prop.begin - subroi_prop.center) / (subroi_prop.sigma*subroi_prop.sigma))*\ + gausFunc(i_t + roi_prop.begin, subroi_prop.center, subroi_prop.sigma, subroi_prop.total_q); + + if (verbose) + std::cout << " tick " << t << " subROI center=" << subroi_prop.center + << " dist=" << sigma_distance << "sigma" + << (sigma_distance <= sigmaWindow ? " [scaled]" : " [pass-through]") << '\n' + << " q_orig=" << q_orig << '\n'; + } + + double delta = q_mod - q_orig; + + // do some sanity checks + if (isnan(q_mod)) { + if (verbose) + std::cout << "WARNING: obtained q_mod = NaN... setting scale to 0" << std::endl; + // for additive correction, no cutoffs are applied + } else if (additiveModification) { + roi_data[i_t] += static_cast(delta); + } else if (q_orig < 0.01) { // check that this is a sane limit + if (verbose) std::cout << "WARNING: obtained q_orig < 0.01 ... setting scale to 1" << std::endl; + } else if (sigma_distance > 9.) { + if (verbose) std::cout << "WARNING: sigma_distance > 9 ... setting scale to 1" << std::endl; + } else { + scale_ratio = q_mod / q_orig; + 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 + << ", scale_ratio=" << scale_ratio << " [multiplicative]" + << ", delta=" << delta << " [additive]" << 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..885a248f0 --- /dev/null +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -0,0 +1,332 @@ +#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 "lardataalg/DetectorInfo/DetectorClocksData.h" +#include "larevt/SpaceCharge/SpaceCharge.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 "sbnobj/ICARUS/TPC/ChannelROI.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? + bool applyXXZAngleScale; // do we scale with X vs XZ angle? + bool applyXdQdXScale; // do we scale with X vs dQ/dX? + bool applyXZAngledQdXScale; // do we scale with XZ angle vs dQ/dX? + bool applyXXWScale; // do we scale with X vs ThXW? + bool additiveModification; // additive (true) vs multiplicative (false) ROI modification + bool useGraph2DInterpolation; // true = Delaunay interp (default), false = nearest bin center + double readoutWindowTicks; // how many ticks are in the readout window? + double tickOffset; // do we want an offset in the ticks? + + 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 + + std::vector graph2Ds_Charge_XXZAngle; // the graphs for the charge correction in X vs XZ angle + std::vector graph2Ds_Sigma_XXZAngle; // the graphs for the width correction in X vs XZ angle + std::vector graph2Ds_Charge_XdQdX; // the graphs for charge correction in X vs dQ/dX + std::vector graph2Ds_Sigma_XdQdX; // the graphs for width correction in X vs dQ/dX + std::vector graph2Ds_Charge_XZAngledQdX; // the graphs for charge correction in XZ angle vs dQ/dX + std::vector graph2Ds_Sigma_XZAngledQdX; // the graphs for width correction in XZ angle vs dQ/dX + std::vector graph2Ds_Charge_XXW; // the graphs for the charge correction in XXW + std::vector graph2Ds_Sigma_XXW; // the graphs for the width correction in XXW + + // lets try making a constructor here + // assume we can get a geometry service, a detector clcok, and a detector properties + // 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 bool& arg_ApplyXXZAngleScale = false, + const bool& arg_ApplyXdQdXScale = false, + const bool& arg_ApplyXZAngledQdXScale = false, + const bool& arg_ApplyXXWScale = true, + const double& arg_TickOffset = 0) + : geometry(geom), + wireReadout(wireRead), + detPropData(detProp), + applyChannelScale(arg_ApplyChannelScale), + applyXScale(arg_ApplyXScale), + applyYZScale(arg_ApplyYZScale), + applyXZAngleScale(arg_ApplyXZAngleScale), + applyYZAngleScale(arg_ApplyYZAngleScale), + applydEdXScale(arg_ApplydEdXScale), + applyXXZAngleScale(arg_ApplyXXZAngleScale), + applyXdQdXScale(arg_ApplyXdQdXScale), + applyXZAngledQdXScale(arg_ApplyXZAngledQdXScale), + applyXXWScale(arg_ApplyXXWScale), + additiveModification(false), + useGraph2DInterpolation(true), + readoutWindowTicks(detProp.ReadOutWindowSize()), // the default A2795 (ICARUS TPC readout board) readout window is 4096 samples + tickOffset(arg_TickOffset) // tick offset is for MC truth, default to zero and set only as necessary + { + } + + // 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; // On Plane0!! + float tick_max; // On Plane0!! + float y; + float z; + double dxdr; + double dydr; + double dzdr; + double dqdr; + double dedr; + ScaleValues_t scales_avg[3]; + } TruthProperties_t; + + // A single SimChannel ionization deposit (sim::IDE) tagged with its readout + // coordinates. This is the SimChannel analogue of a sim::SimEnergyDeposit used by + // the truth-matching path; unlike an edep, a sim::IDE carries only a point position + // (x,y,z) -- direction/step come from the matched simb::MCParticle trajectory. + typedef struct MatchedIDE + { + raw::ChannelID_t channel; // readout channel the charge landed on + geo::WireID wire; // wire closest to the deposit + double tick; // readout tick (TDC converted to tick) + const sim::IDE* ide; // trackID, numElectrons, energy, x, y, z + } MatchedIDE_t; + + // internal containers + std::map< ROI_Key_t, std::vector > ROIMatchedEdepMap; + std::map< ROI_Key_t, std::vector > ROIMatchedHitMap; + + // SimChannel/IDE truth-matching containers (parallel to the Edep versions). + // fIDEVec is the flat owner; ROIMatchedIDEMap stores indices into it, keyed by ROI. + std::vector fIDEVec; + std::map< ROI_Key_t, std::vector > ROIMatchedIDEMap; + + // Space-charge provider, set by the module (lar::providerFrom()). + // Used to map IDE (at-the-wire) positions to/from the true MCParticle trajectory + // frame. When null or SCE disabled, the mappings are the identity. + const spacecharge::SpaceCharge* fSCE = nullptr; + + // some useful functions + // geometries + // TODO is this the most efficient for new v10 iterators? + 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 sinPlaneAngle = std::sin(planeAngle); + double cosPlaneAngle = std::cos(planeAngle); + + //double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; // don't need to rotate Y for this angle + double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; + + double theta = std::atan2(dxdr, dzdrPlaneRel); + return FoldAngle(theta); + } + + double ThetaYZ_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) + { + double sinPlaneAngle = std::sin(planeAngle); + double cosPlaneAngle = std::cos(planeAngle); + + 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 + + double ThetaXY_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) + { + double sinPlaneAngle = std::sin(planeAngle); + double cosPlaneAngle = std::cos(planeAngle); + + double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; + //double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; // don't need to rotate Z for this angle + + double theta = std::atan2(dxdr, dydrPlaneRel); + return FoldAngle(theta); + } + + double ThetaXW(double dxdr, double dydr, double dzdr, double planeAngle) + { + // planeAngle is the wire angle from +z (PlaneGeo::ThetaZ, numerically equal to + // WireReadoutGeom::WireAngleToVertical). The pitch "gamma" angle is measured from the + // wire-normal, i.e. planeAngle - pi/2 -- the same convention as the canonical LArSoft + // pitch (Calorimetry / TrackCaloSkimmer use WireAngleToVertical - pi/2). Subtracting a + // full pi projected the direction onto the wire itself (a cos<->sin swap), yielding the + // wrong angle; subtract pi/2 so cosG is the projection onto the wire-normal (pitch). + double sin = std::sin(planeAngle - 0.5 * ::util::pi<>()); + double cos = std::cos(planeAngle - 0.5 * ::util::pi<>()); + + double cosG = std::abs(dydr * sin + dzdr * cos); + double theta = std::atan(dxdr / cosG); + return std::abs(theta); + } + + // these are set in the .cc file + ROIProperties_t CalcROIProperties(recob::Wire const&, size_t const&); + ROIProperties_t CalcROIProperties(recob::ChannelROI const&, size_t const&); + + std::vector> GetTargetROIs(sim::SimEnergyDeposit const&, double offset); + std::vector> GetHitTargetROIs(recob::Hit const&); + + void FillROIMatchedEdepMap(std::vector const&, std::vector const&, double offset); + void FillROIMatchedHitMap(std::vector const&, std::vector const&); + + // SimChannel/IDE versions of FillROIMatchedEdepMap. The DetectorClocksData is used to + // convert each IDE's TDC to a readout tick (no X->tick projection / tickOffset needed, + // since the SimChannel already carries the readout coordinate). + void FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, detinfo::DetectorClocksData const& clockData, double offset = 0); + void FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& chanROIVec, detinfo::DetectorClocksData const& clockData, double offset = 0); + + void FillROIMatchedEdepMap(std::vector const&, std::vector const&, double offset); + void FillROIMatchedHitMap(std::vector const&, std::vector const&); + + std::vector CalcSubROIProperties(ROIProperties_t const&, std::vector const&); + + std::map> MatchEdepsToSubROIs(std::vector const&, std::vector const&, double offset); + + // SimChannel/IDE analogue of MatchEdepsToSubROIs. Keyed on the IDE's native readout + // tick (no per-plane X->tick projection), otherwise the same center+/-sigma / closest + // sub-ROI assignment as the edep version. + std::map> MatchIDEsToSubROIs(std::vector const&, std::vector const&); + + TruthProperties_t CalcPropertiesFromEdeps(std::vector const&, double offset); + + // SimChannel/IDE analogue of CalcPropertiesFromEdeps. Fills TruthProperties_t from the + // dominant track's sim::IDEs (charge-weighted position/energy) plus the matched + // simb::MCParticle trajectory (direction / pitch / dedr / dqdr). particleMap maps + // abs(trackID) -> MCParticle so the dominant track's trajectory can be looked up. + TruthProperties_t CalcPropertiesFromIDEs(std::vector const&, std::map const& particleMap); + + // Space-charge coordinate maps (mirroring TrackCaloSkimmer). Identity if fSCE is null + // or SCE disabled. + geo::Point_t WireToTrajectoryPosition(const geo::Point_t& loc, const geo::TPCID& tpc) const; // at-the-wire -> true trajectory frame + geo::Point_t TrajectoryToWirePosition(const geo::Point_t& loc, const geo::Vector_t& driftdir) const; // true -> at-the-wire frame + + ScaleValues_t GetScaleValues(TruthProperties_t const&, ROIProperties_t const&); + ScaleValues_t GetChannelScaleValues(TruthProperties_t const&, raw::ChannelID_t const&); + ScaleValues_t GetViewScaleValues(TruthProperties_t const&, geo::View_t const&); + + void ModifyROI(std::vector &, + ROIProperties_t const &, + std::vector const&, + std::map const&, + double sigmaWindow = std::numeric_limits::infinity()); + }; // end class +} // end namespace + +#endif // sbncode_WireMod_Utility_WireModUtility_hh