diff --git a/sbncode/CMakeLists.txt b/sbncode/CMakeLists.txt index 9b7cfc5ed..daed2590a 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(HitFinder) # Supera # diff --git a/sbncode/HitFinder/CMakeLists.txt b/sbncode/HitFinder/CMakeLists.txt new file mode 100644 index 000000000..1d4870cef --- /dev/null +++ b/sbncode/HitFinder/CMakeLists.txt @@ -0,0 +1,58 @@ +add_subdirectory(HitFinderUtilities) + +art_make_library( + LIBRARIES + lardataobj::RawData + lardataobj::RecoBase + lardata::Utilities + fhiclcpp::fhiclcpp + cetlib::cetlib +) +set( MODULE_LIBRARIES + larcorealg::Geometry + larcorealg::Geometry + larcore::Geometry_Geometry_service + lardata::Utilities + larevt::Filters + lardataobj::RawData + larevt::CalibrationDBI_IOVData + larevt::CalibrationDBI_Providers + lardataobj::RecoBase + lardata::ArtDataHelper + larreco::RecoAlg + sbnobj::ICARUS_TPC + sbnobj::Common_Utilities + sbncode::HitFinder_HitFinderUtilities + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art_root_io::tfile_support ROOT::Core + art_root_io::TFileService_service + art::Persistency_Common + art::Persistency_Provenance + art::Utilities + canvas::canvas + messagefacility::MF_MessageLogger + messagefacility::headers + fhiclcpp::fhiclcpp + cetlib::cetlib + ROOT::Geom + ROOT::XMLIO + ROOT::Gdml + ROOT::FFTW + + ) +cet_build_plugin(GaussHitFinderSBN art::module + LIBRARIES ${MODULE_LIBRARIES} + larreco::HitFinder + larreco::CandidateHitFinderTool + larreco::PeakFitterTool + ) +cet_build_plugin(ChannelROIToWire art::module LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(WireToChannelROI art::module LIBRARIES ${MODULE_LIBRARIES}) + + +install_headers() +install_fhicl() +install_source() + diff --git a/sbncode/HitFinder/ChannelROIToWire_module.cc b/sbncode/HitFinder/ChannelROIToWire_module.cc new file mode 100644 index 000000000..243d5a744 --- /dev/null +++ b/sbncode/HitFinder/ChannelROIToWire_module.cc @@ -0,0 +1,166 @@ +//////////////////////////////////////////////////////////////////////// +// +// ChannelROIToWire class - Converts from ChannelROI to Wire data +// +// usher@slac.stanford.edu +// +//////////////////////////////////////////////////////////////////////// + +// C/C++ standard libraries +#include +#include +#include // std::pair<> +#include // std::unique_ptr<> + +// framework libraries +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "canvas/Utilities/Exception.h" +#include "canvas/Utilities/InputTag.h" + +// LArSoft libraries +#include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t +#include "larcore/Geometry/WireReadout.h" +#include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom() +#include "larcorealg/CoreUtils/zip.h" +#include "lardataobj/RecoBase/Wire.h" +#include "lardata/ArtDataHelper/WireCreator.h" + +#include "sbnobj/ICARUS/TPC/ChannelROI.h" + +///creation of calibrated signals on wires +namespace caldata { + +class ChannelROIToWire : public art::EDProducer +{ +public: +// create calibrated signals on wires. this class runs +// an fft to remove the electronics shaping. + explicit ChannelROIToWire(fhicl::ParameterSet const& pset); + void produce(art::Event& evt); + void beginJob(); + void endJob(); + void reconfigure(fhicl::ParameterSet const& p); +private: + + std::vector fWireModuleLabelVec; ///< vector of modules that made digits + std::vector fOutInstanceLabelVec; ///< The output instance labels to apply + bool fDiagnosticOutput; ///< secret diagnostics flag + size_t fEventCount; ///< count of event processed + + const geo::WireReadoutGeom* fChannelMapAlg = &art::ServiceHandle()->Get(); + +}; // class ChannelROIToWire + +DEFINE_ART_MODULE(ChannelROIToWire) + +//------------------------------------------------- +ChannelROIToWire::ChannelROIToWire(fhicl::ParameterSet const& pset) : EDProducer{pset} +{ + this->reconfigure(pset); + + for(const auto& wireLabel : fOutInstanceLabelVec) + { + produces>(wireLabel); + } +} + +////////////////////////////////////////////////////// +void ChannelROIToWire::reconfigure(fhicl::ParameterSet const& pset) +{ + // Recover the parameters + fWireModuleLabelVec = pset.get>("WireModuleLabelVec", std::vector()={"decon1droi"}); + fOutInstanceLabelVec = pset.get> ("OutInstanceLabelVec", {"PHYSCRATEDATA"}); + fDiagnosticOutput = pset.get< bool >("DiagnosticOutput", false); + + if (fWireModuleLabelVec.size() != fOutInstanceLabelVec.size()) + { + throw art::Exception(art::errors::Configuration) << " Configured " << fOutInstanceLabelVec.size() + << " instance names (`OutInstanceLabelVec`) for " << fWireModuleLabelVec.size() + << " input products (`WireModuleLabelVec`)\n"; + } + + return; +} + +//------------------------------------------------- +void ChannelROIToWire::beginJob() +{ + fEventCount = 0; +} // beginJob + +////////////////////////////////////////////////////// +void ChannelROIToWire::endJob() +{ +} + +////////////////////////////////////////////////////// +void ChannelROIToWire::produce(art::Event& evt) +{ + // We need to loop through the list of Wire data we have been given + // This construct from Gianluca Petrillo who invented it and should be given all credit for it! + for(auto const& [channelLabel, instanceName] : util::zip(fWireModuleLabelVec, fOutInstanceLabelVec)) + { + // make a collection of Wires + std::unique_ptr> wireCol = std::make_unique>(); + + mf::LogInfo("ChannelROIToWire") << "ChannelROIToWire, looking for ChannelROI data at " << channelLabel.encode(); + + // Read in the collection of full length deconvolved waveforms + const std::vector& channelVec = evt.getProduct>(channelLabel); + + mf::LogInfo("ChannelROIToWire") << "--> Recovered ChannelROI data, size: " << channelVec.size(); + + if (!channelVec.empty()) + { + // Reserve the room for the output + wireCol->reserve(channelVec.size()); + + // Loop through the input ChannelROI collection + for(const auto& channelROI : channelVec) + { + // Recover the channel and the view + raw::ChannelID_t channel = channelROI.Channel(); + geo::View_t view = fChannelMapAlg->View(channel); + float ADCScaleFactor = channelROI.ADCScaleFactor(); + + // 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()) + { + 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()); + } + + mf::LogInfo("ChannelROIToWire") << "--> Outputting Wire data, size: " << wireCol->size() << " with instance name: " << instanceName; + + // Time to stroe everything + if(wireCol->empty()) mf::LogWarning("ChannelROIToWire") << "No wires made for this event."; + } + + evt.put(std::move(wireCol), instanceName); + } + + fEventCount++; + + return; +} // produce + +} // end namespace caldata diff --git a/sbncode/HitFinder/GaussHitFinderSBN_module.cc b/sbncode/HitFinder/GaussHitFinderSBN_module.cc new file mode 100644 index 000000000..c0a91cc91 --- /dev/null +++ b/sbncode/HitFinder/GaussHitFinderSBN_module.cc @@ -0,0 +1,669 @@ +//////////////////////////////////////////////////////////////////////// +// +// GaussHitFinder class +// +// jaasaadi@syr.edu +// +// This algorithm is designed to find hits on wires after deconvolution. +// ----------------------------------- +// This algorithm is based on the FFTHitFinder written by Brian Page, +// Michigan State University, for the ArgoNeuT experiment. +// +// +// The algorithm walks along the wire and looks for pulses above threshold +// The algorithm then attempts to fit n-gaussians to these pulses where n +// is set by the number of peaks found in the pulse +// If the Chi2/NDF returned is "bad" it attempts to fit n+1 gaussians to +// the pulse. If this is a better fit it then uses the parameters of the +// Gaussian fit to characterize the "hit" object +// +// To use this simply include the following in your producers: +// gaushit: @local::microboone_GaussHitFinderSBN +// gaushit: @local::argoneut_GaussHitFinderSBN +//////////////////////////////////////////////////////////////////////// + +// C/C++ standard library +#include // std::accumulate() +#include +#include // std::unique_ptr() +#include +#include // std::move() + +// Framework includes +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Core/SharedProducer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art/Utilities/Globals.h" +#include "art/Utilities/make_tool.h" +#include "art_root_io/TFileService.h" +#include "canvas/Persistency/Common/FindOneP.h" +#include "fhiclcpp/ParameterSet.h" + +// LArSoft Includes +#include "larcore/Geometry/WireReadout.h" +#include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t +#include "larreco/HitFinder/HitFilterAlg.h" +#include "lardataobj/RecoBase/Hit.h" + +#include "sbnobj/ICARUS/TPC/ChannelROI.h" +#include "sbncode/HitFinder/HitFinderUtilities/HitCreator.h" + +#include "larreco/HitFinder/HitFinderTools/ICandidateHitFinder.h" +#include "larreco/HitFinder/HitFinderTools/IPeakFitter.h" + +// ROOT Includes +#include "TH1F.h" +#include "TMath.h" + +#include "tbb/concurrent_vector.h" +#include "tbb/parallel_for.h" + +namespace hit { + class GaussHitFinderSBN : public art::SharedProducer { + public: + explicit GaussHitFinderSBN(fhicl::ParameterSet const& pset, art::ProcessingFrame const&); + + private: + void produce(art::Event& evt, art::ProcessingFrame const&) override; + void beginJob(art::ProcessingFrame const&) override; + + std::vector FillOutHitParameterVector(const std::vector& input); + + const bool fFilterHits; + const bool fFillHists; + + const std::string fCalDataModuleLabel; + const std::string fAllHitsInstanceName; + + const std::vector fLongMaxHitsVec; /// fLongPulseWidthVec; /// + fAreaNormsVec; /// fPulseHeightCuts; + const std::vector fPulseWidthCuts; + const std::vector fPulseRatioCuts; + + std::atomic fEventCount{0}; + + //only Standard and Morphological implementation is threadsafe. + std::vector> + fHitFinderToolVec; ///< For finding candidate hits + // only Marqdt implementation is threadsafe. + std::unique_ptr fPeakFitterTool; ///< Perform fit to candidate peaks + //HitFilterAlg implementation is threadsafe. + std::unique_ptr fHitFilterAlg; ///< algorithm used to filter out noise hits + + //only used when fFillHists is true and in single threaded mode. + TH1F* fFirstChi2; + TH1F* fChi2; + + }; // class GaussHitFinderSBN + + //------------------------------------------------- + //------------------------------------------------- + GaussHitFinderSBN::GaussHitFinderSBN(fhicl::ParameterSet const& pset, art::ProcessingFrame const&) + : SharedProducer{pset} + , fFilterHits(pset.get("FilterHits", false)) + , fFillHists(pset.get("FillHists", false)) + , fCalDataModuleLabel(pset.get("CalDataModuleLabel")) + , fAllHitsInstanceName(pset.get("AllHitsInstanceName", "")) + , fLongMaxHitsVec(pset.get>("LongMaxHits", std::vector() = {25, 25, 25})) + , fLongPulseWidthVec( + pset.get>("LongPulseWidth", std::vector() = {16, 16, 16})) + , fMaxMultiHit(pset.get("MaxMultiHit")) + , fAreaMethod(pset.get("AreaMethod")) + , fAreaNormsVec(FillOutHitParameterVector(pset.get>("AreaNorms"))) + , fChi2NDF(pset.get("Chi2NDF")) + , fPulseHeightCuts( + pset.get>("PulseHeightCuts", std::vector() = {3.0, 3.0, 3.0})) + , fPulseWidthCuts( + pset.get>("PulseWidthCuts", std::vector() = {2.0, 1.5, 1.0})) + , fPulseRatioCuts( + pset.get>("PulseRatioCuts", std::vector() = {0.35, 0.40, 0.20})) + { + if (fFillHists && art::Globals::instance()->nthreads() > 1u) { + throw art::Exception(art::errors::Configuration) + << "Cannot fill histograms when multiple threads configured, please set fFillHists to " + "false or change number of threads to 1\n"; + } + async(); + if (fFilterHits) { + fHitFilterAlg = std::make_unique(pset.get("HitFilterAlg")); + } + + // recover the tool to do the candidate hit finding + // Recover the vector of fhicl parameters for the ROI tools + const fhicl::ParameterSet& hitFinderTools = pset.get("HitFinderToolVec"); + + fHitFinderToolVec.resize(hitFinderTools.get_pset_names().size()); + + for (const std::string& hitFinderTool : hitFinderTools.get_pset_names()) { + const fhicl::ParameterSet& hitFinderToolParamSet = + hitFinderTools.get(hitFinderTool); + size_t planeIdx = hitFinderToolParamSet.get("Plane"); + + fHitFinderToolVec.at(planeIdx) = + art::make_tool(hitFinderToolParamSet); + } + + // Recover the peak fitting tool + fPeakFitterTool = + art::make_tool(pset.get("PeakFitter")); + + // let HitCollectionCreator declare that we are going to produce + // hits and associations with wires and raw digits + // We want the option to output two hit collections, one filtered + // and one with all hits. The key to doing this will be a non-null + // instance name for the second collection + // (with no particular product label) + sbn::HitCollectionCreator::declare_products( + producesCollector(), fAllHitsInstanceName, true, false); //fMakeRawDigitAssns); + + // and now the filtered hits... + if (fAllHitsInstanceName != "") + sbn::HitCollectionCreator::declare_products( + producesCollector(), "", true, false); //fMakeRawDigitAssns); + + return; + } // GaussHitFinderSBN::GaussHitFinderSBN() + + //------------------------------------------------- + //------------------------------------------------- + std::vector GaussHitFinderSBN::FillOutHitParameterVector(const std::vector& input) + { + if (input.size() == 0) + throw std::runtime_error( + "GaussHitFinderSBN::FillOutHitParameterVector ERROR! Input config vector has zero size."); + + std::vector output; + auto const& wireReadoutGeom = art::ServiceHandle()->Get(); + const unsigned int N_PLANES = wireReadoutGeom.Nplanes(); + + if (input.size() == 1) + output.resize(N_PLANES, input[0]); + else if (input.size() == N_PLANES) + output = input; + else + throw std::runtime_error("GaussHitFinderSBN::FillOutHitParameterVector ERROR! Input config " + "vector size !=1 and !=N_PLANES."); + return output; + } + + //------------------------------------------------- + //------------------------------------------------- + void GaussHitFinderSBN::beginJob(art::ProcessingFrame const&) + { + // get access to the TFile service + art::ServiceHandle tfs; + + // ====================================== + // === Hit Information for Histograms === + if (fFillHists) { + fFirstChi2 = tfs->make("fFirstChi2", "#chi^{2}", 10000, 0, 5000); + fChi2 = tfs->make("fChi2", "#chi^{2}", 10000, 0, 5000); + } + } + + // This algorithm uses the fact that deconvolved signals are very smooth + // and looks for hits as areas between local minima that have signal above + // threshold. + //------------------------------------------------- + void GaussHitFinderSBN::produce(art::Event& evt, art::ProcessingFrame const&) + { + unsigned int count = fEventCount.fetch_add(1); + //================================================================================================== + + TH1::AddDirectory(kFALSE); + + auto const& wireReadoutGeom = art::ServiceHandle()->Get(); + + // ############################################### + // ### Making a ptr vector to put on the event ### + // ############################################### + // this contains the hit collection + // and its associations to wires and raw digits + sbn::HitCollectionCreator allHitCol(evt, fAllHitsInstanceName, true, false); + + // Handle the filtered hits collection... + sbn::HitCollectionCreator hcol(evt, "", true, false); + sbn::HitCollectionCreator* filteredHitCol = 0; + + if (fFilterHits) filteredHitCol = &hcol; + + //store in a thread safe way + struct hitstruct { + recob::Hit hit_tbb; + art::Ptr channelROI_tbb; + }; + + tbb::concurrent_vector hitstruct_vec; + tbb::concurrent_vector filthitstruct_vec; + + // ########################################## + // ### Reading in the Wire List object(s) ### + // ########################################## + art::Handle> channelROIVecHandle; + evt.getByLabel(fCalDataModuleLabel, channelROIVecHandle); + + //################################################# + //### Set the charge determination method ### + //### Default is to compute the normalized area ### + //################################################# + std::function chargeFunc = + [](double /* peakMean */, + double peakAmp, + double peakWidth, + double areaNorm, + int /* low */, + int /* hi */) { return std::sqrt(2 * TMath::Pi()) * peakAmp * peakWidth / areaNorm; }; + + //############################################## + //### Alternative is to integrate over pulse ### + //############################################## + if (fAreaMethod == 0) + chargeFunc = [](double peakMean, + double peakAmp, + double peakWidth, + double /* areaNorm */, + int low, + int hi) { + double charge(0); + for (int sigPos = low; sigPos < hi; sigPos++) + charge += peakAmp * TMath::Gaus(sigPos, peakMean, peakWidth); + return charge; + }; + + //############################## + //### Looping over the wires ### + //############################## + tbb::parallel_for( + static_cast(0), + channelROIVecHandle->size(), + [&](size_t& channelROIIter) { + // #################################### + // ### Getting this particular channelROI ### + // #################################### + art::Ptr channelROI(channelROIVecHandle, channelROIIter); + + // --- Setting Channel Number and Signal type --- + + raw::ChannelID_t channel = channelROI->Channel(); + + // get the WireID for this hit + std::vector wids = wireReadoutGeom.ChannelToWire(channel); + // for now, just take the first option returned from ChannelToWire + geo::WireID wid = wids[0]; + // We need to know the plane to look up parameters + geo::PlaneID::PlaneID_t plane = wid.Plane; + + // ---------------------------------------------------------- + // -- Setting the appropriate signal widths and thresholds -- + // -- for the right plane. -- + // ---------------------------------------------------------- + + // ################################################# + // ### Set up to loop over ROI's for this wire ### + // ################################################# + const recob::ChannelROI::RegionsOfInterest_f signalROI = channelROI->SignalROIF(); + + tbb::parallel_for( + static_cast(0), + signalROI.n_ranges(), + [&](size_t& rangeIter) { + const auto& range = signalROI.range(rangeIter); + // ROI start time + raw::TDCtick_t roiFirstBinTick = range.begin_index(); + + // ########################################################### + // ### Scan the waveform and find candidate peaks + merge ### + // ########################################################### + + reco_tool::ICandidateHitFinder::HitCandidateVec hitCandidateVec; + reco_tool::ICandidateHitFinder::MergeHitCandidateVec mergedCandidateHitVec; + + fHitFinderToolVec.at(plane)->findHitCandidates( + range, 0, channel, count, hitCandidateVec); + fHitFinderToolVec.at(plane)->MergeHitCandidates( + range, hitCandidateVec, mergedCandidateHitVec); + + // ####################################################### + // ### Lets loop over the pulses we found on this wire ### + // ####################################################### + + for (auto& mergedCands : mergedCandidateHitVec) { + int startT = mergedCands.front().startTick; + int endT = mergedCands.back().stopTick; + + // ### Putting in a protection in case things went wrong ### + // ### In the end, this primarily catches the case where ### + // ### a fake pulse is at the start of the ROI ### + if (endT - startT < 5) continue; + + // ####################################################### + // ### Clearing the parameter vector for the new Pulse ### + // ####################################################### + + // === Setting The Number Of Gaussians to try === + int nGausForFit = mergedCands.size(); + + // ################################################## + // ### Calling the function for fitting Gaussians ### + // ################################################## + double chi2PerNDF(0.); + int NDF(1); + /*stand alone + reco_tool::IPeakFitter::PeakParamsVec peakParamsVec(nGausForFit); + */ + reco_tool::IPeakFitter::PeakParamsVec peakParamsVec; + + // ####################################################### + // ### If # requested Gaussians is too large then punt ### + // ####################################################### + if (mergedCands.size() <= fMaxMultiHit) { + fPeakFitterTool->findPeakParameters( + range.data(), mergedCands, peakParamsVec, chi2PerNDF, NDF); + + // If the chi2 is infinite then there is a real problem so we bail + if (!(chi2PerNDF < std::numeric_limits::infinity())) { + chi2PerNDF = 2. * fChi2NDF; + NDF = 2; + } + + if (fFillHists) fFirstChi2->Fill(chi2PerNDF); + } + + // ####################################################### + // ### If too large then force alternate solution ### + // ### - Make n hits from pulse train where n will ### + // ### depend on the fhicl parameter fLongPulseWidth ### + // ### Also do this if chi^2 is too large ### + // ####################################################### + if (mergedCands.size() > fMaxMultiHit || nGausForFit * chi2PerNDF > fChi2NDF) { + int longPulseWidth = fLongPulseWidthVec.at(plane); + int nHitsThisPulse = (endT - startT) / longPulseWidth; + + if (nHitsThisPulse > fLongMaxHitsVec.at(plane)) { + nHitsThisPulse = fLongMaxHitsVec.at(plane); + longPulseWidth = (endT - startT) / nHitsThisPulse; + } + + if (nHitsThisPulse * longPulseWidth < endT - startT) nHitsThisPulse++; + + int firstTick = startT; + int lastTick = std::min(firstTick + longPulseWidth, endT); + + peakParamsVec.clear(); + nGausForFit = nHitsThisPulse; + NDF = 1.; + chi2PerNDF = chi2PerNDF > fChi2NDF ? chi2PerNDF : -1.; + + for (int hitIdx = 0; hitIdx < nHitsThisPulse; hitIdx++) { + // This hit parameters + double ROIsumADC = + std::accumulate(range.begin() + firstTick, range.begin() + lastTick, 0.); + double peakSigma = (lastTick - firstTick) / 3.; // Set the width... + double peakAmp = 0.3989 * ROIsumADC / peakSigma; // Use gaussian formulation + double peakMean = (firstTick + lastTick) / 2.; + + // Store hit params + reco_tool::IPeakFitter::PeakFitParams_t peakParams; + + peakParams.peakCenter = peakMean; + peakParams.peakCenterError = 0.1 * peakMean; + peakParams.peakSigma = peakSigma; + peakParams.peakSigmaError = 0.1 * peakSigma; + peakParams.peakAmplitude = peakAmp; + peakParams.peakAmplitudeError = 0.1 * peakAmp; + + peakParamsVec.push_back(peakParams); + + // set for next loop + firstTick = lastTick; + lastTick = std::min(lastTick + longPulseWidth, endT); + } + } + + // ####################################################### + // ### Loop through returned peaks and make recob hits ### + // ####################################################### + + int numHits(0); + + // Make a container for what will be the filtered collection + std::vector filteredHitVec; + + float nextpeak(0); + float prevpeak(0); + float nextpeakSig(0); + float prevpeakSig(0); + float nsigmaADC(2.0); + float newright(0); + float newleft(0); + for (const auto& peakParams : peakParamsVec) { + // Extract values for this hit + float peakAmp = peakParams.peakAmplitude; + float peakMean = peakParams.peakCenter; + float peakWidth = peakParams.peakSigma; + + //std::cout<<" ans hits "< 0) { + prevpeak = (peakParamsVec.at(numHits - 1)).peakCenter; + prevpeakSig = (peakParamsVec.at(numHits - 1)).peakSigma; + //std::cout<<" ans size "<::const_iterator sumStartItr = range.begin() + startT; + std::vector::const_iterator sumEndItr = range.begin() + endT; + + //### limits for the sum of the Hit based on the gaussian peak and sigma + std::vector::const_iterator HitsumStartItr = + range.begin() + peakMean - nsigmaADC * peakWidth; + std::vector::const_iterator HitsumEndItr = + range.begin() + peakMean + nsigmaADC * peakWidth; + + if (nGausForFit > 1) { + if (numHits > 0) { + if ((peakMean - nsigmaADC * peakWidth) < (prevpeak + nsigmaADC * prevpeakSig)) { + float difPeak = peakMean - prevpeak; + float weightpeak = prevpeakSig / (prevpeakSig + peakWidth); + HitsumStartItr = range.begin() + prevpeak + difPeak * weightpeak; + newleft = prevpeak + difPeak * weightpeak; + } + } + + if (numHits < nGausForFit - 1) { + if ((peakMean + nsigmaADC * peakWidth) > (nextpeak - nsigmaADC * nextpeakSig)) { + float difPeak = nextpeak - peakMean; + float weightpeak = peakWidth / (nextpeakSig + peakWidth); + HitsumEndItr = range.begin() + peakMean + difPeak * weightpeak; + newright = peakMean + difPeak * weightpeak; + } + } + } + + //protection to avoid negative ranges + if (newright - newleft < 0) continue; + if (HitsumStartItr > HitsumEndItr) continue; + + //avoid ranges out of ROI if it happens + if (HitsumStartItr < sumStartItr) HitsumStartItr = sumStartItr; + + if (HitsumEndItr > sumEndItr) HitsumEndItr = sumEndItr; + + // ### Sum of ADC counts + float ROIsumADC = std::accumulate(sumStartItr, sumEndItr, 0.); + float HitsumADC = std::accumulate(HitsumStartItr, HitsumEndItr, 0.); + float peakTime = peakMean + roiFirstBinTick; + float chi2ByNDF = chi2PerNDF; // convert to float + + // ok, now create the hit + sbn::HitCreator hitcreator( + *channelROI, // channelROI reference + wid, // wire ID + startT + roiFirstBinTick, // start_tick TODO check + endT + roiFirstBinTick, // end_tick TODO check + peakWidth, // rms + peakTime, // peak_time + peakMeanErr, // sigma_peak_time + peakAmp, // peak_amplitude + peakAmpErr, // sigma_peak_amplitude + charge, // hit_integral + chargeErr, // hit_sigma_integral + ROIsumADC, // summedADC of the ROI + HitsumADC, // summedADC of the Hit + nGausForFit, // multiplicity + numHits, // local_index TODO check that the order is correct + chi2ByNDF, // goodness_of_fit + NDF // dof + ); + + if (filteredHitCol) filteredHitVec.push_back(hitcreator.copy()); + + const recob::Hit hit(hitcreator.move()); + + // This loop will store ALL hits + hitstruct tmp{std::move(hit), channelROI}; + hitstruct_vec.push_back(std::move(tmp)); + + numHits++; + } // <---End loop over gaussians + + // Should we filter hits? + if (filteredHitCol && !filteredHitVec.empty()) { + // ####################################################################### + // Is all this sorting really necessary? Would it be faster to just loop + // through the hits and perform simple cuts on amplitude and width on a + // hit-by-hit basis, either here in the module (using fPulseHeightCuts and + // fPulseWidthCuts) or in HitFilterAlg? + // ####################################################################### + + // Sort in ascending peak height + std::sort(filteredHitVec.begin(), + filteredHitVec.end(), + [](const auto& left, const auto& right) { + return left.PeakAmplitude() > right.PeakAmplitude(); + }); + + // Reject if the first hit fails the PH/wid cuts + if (filteredHitVec.front().PeakAmplitude() < fPulseHeightCuts.at(plane) || + filteredHitVec.front().RMS() < fPulseWidthCuts.at(plane)) + filteredHitVec.clear(); + + // Now check other hits in the snippet + if (filteredHitVec.size() > 1) { + // The largest pulse height will now be at the front... + float largestPH = filteredHitVec.front().PeakAmplitude(); + + // Find where the pulse heights drop below threshold + float threshold(fPulseRatioCuts.at(plane)); + + std::vector::iterator smallHitItr = + std::find_if(filteredHitVec.begin(), + filteredHitVec.end(), + [largestPH, threshold](const auto& hit) { + return hit.PeakAmplitude() < 8. && + hit.PeakAmplitude() / largestPH < threshold; + }); + + // Shrink to fit + if (smallHitItr != filteredHitVec.end()) + filteredHitVec.resize(std::distance(filteredHitVec.begin(), smallHitItr)); + + // Resort in time order + std::sort(filteredHitVec.begin(), + filteredHitVec.end(), + [](const auto& left, const auto& right) { + return left.PeakTime() < right.PeakTime(); + }); + } + + // Copy the hits we want to keep to the filtered hit collection +// for (const auto& filteredHit : filteredHitVec) +// if (!fHitFilterAlg || fHitFilterAlg->IsGoodHit(filteredHit)) { +// hitstruct tmp{std::move(filteredHit), channelROI}; +// filthitstruct_vec.push_back(std::move(tmp)); +// } +// + if (fFillHists) fChi2->Fill(chi2PerNDF); + } + } //<---End loop over merged candidate hits + } //<---End looping over ROI's + ); //end tbb parallel for + } //<---End looping over all the wires + ); //end tbb parallel for + + for (size_t i = 0; i < hitstruct_vec.size(); i++) { + allHitCol.emplace_back(hitstruct_vec[i].hit_tbb, hitstruct_vec[i].channelROI_tbb); + } + + for (size_t j = 0; j < filthitstruct_vec.size(); j++) { + filteredHitCol->emplace_back(filthitstruct_vec[j].hit_tbb, filthitstruct_vec[j].channelROI_tbb); + } + + //================================================================================================== + // End of the event -- move the hit collection and the associations into the event + + if (filteredHitCol) { + + // If we filtered hits but no instance name was + // specified for the "all hits" collection, then + // only save the filtered hits to the event + if (fAllHitsInstanceName == "") { + filteredHitCol->put_into(evt); + + // otherwise, save both + } + else { + filteredHitCol->put_into(evt); + allHitCol.put_into(evt); + } + } + else { + allHitCol.put_into(evt); + } + } // End of produce() + + DEFINE_ART_MODULE(GaussHitFinderSBN) + +} // end of hit namespace diff --git a/sbncode/HitFinder/HitFinderUtilities/CMakeLists.txt b/sbncode/HitFinder/HitFinderUtilities/CMakeLists.txt new file mode 100644 index 000000000..2e9266b58 --- /dev/null +++ b/sbncode/HitFinder/HitFinderUtilities/CMakeLists.txt @@ -0,0 +1,21 @@ +cet_make_library( + SOURCE + HitCreator.cxx + LIBRARIES + cetlib_except::cetlib_except + messagefacility::MF_MessageLogger + sbnobj::ICARUS_TPC + larcore::Geometry_Geometry_service + larcorealg::Geometry + lardata::ArtDataHelper + lardataobj::AnalysisBase + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art::Persistency_Common + art::Persistency_Provenance + art::Utilities + ) + +install_headers() +install_source() diff --git a/sbncode/HitFinder/HitFinderUtilities/HitCreator.cxx b/sbncode/HitFinder/HitFinderUtilities/HitCreator.cxx new file mode 100644 index 000000000..76f9833ae --- /dev/null +++ b/sbncode/HitFinder/HitFinderUtilities/HitCreator.cxx @@ -0,0 +1,638 @@ +/** **************************************************************************** + * @file HitCreator.cxx + * @brief Helper functions to create a hit from ChannelROIs - implementation file + * @date April 20, 2025 + * @author usher@slac.stanford.edu + * @see Hit.h HitCreator.h + * + * ****************************************************************************/ + +// declaration header +#include "sbncode/HitFinder/HitFinderUtilities/HitCreator.h" + +// C/C++ standard library +#include // std::accumulate(), std::max() +#include +#include // std::numeric_limits<> +#include // std::move() + +// art libraries +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "canvas/Persistency/Common/FindOneP.h" +#include "canvas/Utilities/Exception.h" + +// LArSoft libraries +#include "larcore/Geometry/WireReadout.h" +#include "lardata/Utilities/MakeIndex.h" + +namespace { + + /// Erases the content of an association + template + void ClearAssociations(art::Assns& assns) + { + art::Assns empty; + assns.swap(empty); + } // ClearAssociations() + +} // local namespace + +/// Reconstruction base classes +namespace sbn { + + //**************************************************************************** + //*** HitCreator + //---------------------------------------------------------------------- + HitCreator::HitCreator(raw::RawDigit const& digits, + geo::WireID const& wireID, + raw::TDCtick_t start_tick, + raw::TDCtick_t end_tick, + float rms, + float peak_time, + float sigma_peak_time, + float peak_amplitude, + float sigma_peak_amplitude, + float hit_integral, + float hit_sigma_integral, + float ROIsummedADC, + float HitsummedADC, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof) + : hit(digits.Channel(), + start_tick, + end_tick, + peak_time, + sigma_peak_time, + rms, + peak_amplitude, + sigma_peak_amplitude, + ROIsummedADC, + HitsummedADC, + hit_integral, + hit_sigma_integral, + multiplicity, + local_index, + goodness_of_fit, + dof, + art::ServiceHandle()->Get().View(digits.Channel()), + art::ServiceHandle()->Get().SignalType(digits.Channel()), + wireID) + {} // HitCreator::HitCreator(RawDigit) + + //---------------------------------------------------------------------- + HitCreator::HitCreator(recob::ChannelROI const& wire, + geo::WireID const& wireID, + raw::TDCtick_t start_tick, + raw::TDCtick_t end_tick, + float rms, + float peak_time, + float sigma_peak_time, + float peak_amplitude, + float sigma_peak_amplitude, + float hit_integral, + float hit_sigma_integral, + float ROIsummedADC, + float HitsummedADC, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof) + : hit(wire.Channel(), + start_tick, + end_tick, + peak_time, + sigma_peak_time, + rms, + peak_amplitude, + sigma_peak_amplitude, + ROIsummedADC, + HitsummedADC, + hit_integral, + hit_sigma_integral, + multiplicity, + local_index, + goodness_of_fit, + dof, + art::ServiceHandle()->Get().View(wire.Channel()), + art::ServiceHandle()->Get().SignalType(wire.Channel()), + wireID) + {} // HitCreator::HitCreator(Wire) + + //---------------------------------------------------------------------- + HitCreator::HitCreator(recob::ChannelROI const& wire, + geo::WireID const& wireID, + raw::TDCtick_t start_tick, + raw::TDCtick_t end_tick, + float rms, + float peak_time, + float sigma_peak_time, + float peak_amplitude, + float sigma_peak_amplitude, + float hit_integral, + float hit_sigma_integral, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof) + : HitCreator(wire, + wireID, + start_tick, + end_tick, + rms, + peak_time, + sigma_peak_time, + peak_amplitude, + sigma_peak_amplitude, + hit_integral, + hit_sigma_integral, + std::accumulate(wire.SignalROI().begin() + start_tick, + wire.SignalROI().begin() + end_tick, + 0.), // sum of ADC counts between start_tick and end_tic + std::accumulate( + wire.SignalROI().begin() + start_tick, + wire.SignalROI().begin() + end_tick, + 0.), // sum of ADC counts between start_tick and end_tick TO BE MODIFIED IF HIT HAS TO BE CONSIDERED + multiplicity, + local_index, + goodness_of_fit, + dof) + {} // HitCreator::HitCreator(Wire; no summed ADC) + + //---------------------------------------------------------------------- + HitCreator::HitCreator(recob::ChannelROI const& wire, + geo::WireID const& wireID, + float rms, + float peak_time, + float sigma_peak_time, + float peak_amplitude, + float sigma_peak_amplitude, + float hit_integral, + float hit_sigma_integral, + float ROIsummedADC, + float HitsummedADC, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof, + RegionOfInterest_t const& signal) + : HitCreator(wire, + wireID, + signal.begin_index(), + signal.end_index(), + rms, + peak_time, + sigma_peak_time, + peak_amplitude, + sigma_peak_amplitude, + hit_integral, + hit_sigma_integral, + ROIsummedADC, + HitsummedADC, + multiplicity, + local_index, + goodness_of_fit, + dof) + {} // HitCreator::HitCreator(Wire; RoI) + + //---------------------------------------------------------------------- + HitCreator::HitCreator(recob::ChannelROI const& wire, + geo::WireID const& wireID, + float rms, + float peak_time, + float sigma_peak_time, + float peak_amplitude, + float sigma_peak_amplitude, + float hit_integral, + float hit_sigma_integral, + float ROIsummedADC, + float HitsummedADC, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof, + size_t iSignalRoI) + : HitCreator(wire, + wireID, + rms, + peak_time, + sigma_peak_time, + peak_amplitude, + sigma_peak_amplitude, + hit_integral, + hit_sigma_integral, + ROIsummedADC, + HitsummedADC, + multiplicity, + local_index, + goodness_of_fit, + dof, + wire.SignalROI().range(iSignalRoI)) + {} // HitCreator::HitCreator(Wire; RoI index) + + HitCreator::HitCreator(recob::Hit const& from) : hit(from) {} + + // Need to do this way since we don't have friendship for the wireID + HitCreator::HitCreator(recob::Hit const& from, geo::WireID const& wireID) + : hit(from.Channel(), + from.StartTick(), + from.EndTick(), + from.PeakTime(), + from.SigmaPeakTime(), + from.RMS(), + from.PeakAmplitude(), + from.SigmaPeakAmplitude(), + from.ROISummedADC(), + from.HitSummedADC(), + from.Integral(), + from.SigmaIntegral(), + from.Multiplicity(), + from.LocalIndex(), + from.GoodnessOfFit(), + from.DegreesOfFreedom(), + art::ServiceHandle()->Get().View(from.Channel()), + art::ServiceHandle()->Get().SignalType(from.Channel()), + wireID) + { } // HitCreator::HitCreator(new wire ID) + + //**************************************************************************** + //*** HitAndAssociationsWriterBase + //---------------------------------------------------------------------- + HitAndAssociationsWriterBase::HitAndAssociationsWriterBase(art::Event& event, + std::string instance_name, + bool doWireAssns, + bool doRawDigitAssns) + : prod_instance(instance_name) + , hits() + , WireAssns(doWireAssns ? new art::Assns : nullptr) + , RawDigitAssns(doRawDigitAssns ? new art::Assns : nullptr) + , event(&event) + , hitPtrMaker(*(this->event), prod_instance) + {} // HitAndAssociationsWriterBase::HitAndAssociationsWriterBase() + + //------------------------------------------------------------------------------ + void HitAndAssociationsWriterBase::declare_products(art::ProducesCollector& collector, + std::string instance_name /* = "" */, + bool doWireAssns /* = true */, + bool doRawDigitAssns /* = true */ + ) + { + collector.produces>(instance_name); + + // declare the other products we are creating (if any) + if (doWireAssns) { collector.produces>(instance_name); } + if (doRawDigitAssns) { + collector.produces>(instance_name); + } + } // HitAndAssociationsWriterBase::declare_products() + + //------------------------------------------------------------------------------ + void HitAndAssociationsWriterBase::put_into() + { + assert(event); + if (hits) event->put(std::move(hits), prod_instance); + if (WireAssns) event->put(std::move(WireAssns), prod_instance); + if (RawDigitAssns) event->put(std::move(RawDigitAssns), prod_instance); + } // HitAndAssociationsWriterBase::put_into() + + //**************************************************************************** + //*** HitCollectionCreator + //---------------------------------------------------------------------- + HitCollectionCreator::HitCollectionCreator(art::Event& event, + std::string instance_name /* = "" */, + bool doWireAssns /* = true */, + bool doRawDigitAssns /* = true */ + ) + : HitAndAssociationsWriterBase(event, instance_name, doWireAssns, doRawDigitAssns) + { + hits.reset(new std::vector); + } // HitCollectionCreator::HitCollectionCreator() + + //---------------------------------------------------------------------- + void HitCollectionCreator::emplace_back(recob::Hit&& hit, + art::Ptr const& wire, + art::Ptr const& digits) + { + + // add the hit to the collection + hits->emplace_back(std::move(hit)); + + CreateAssociationsToLastHit(wire, digits); + } // HitCollectionCreator::emplace_back(Hit&&) + + //---------------------------------------------------------------------- + void HitCollectionCreator::emplace_back(recob::Hit const& hit, + art::Ptr const& wire, + art::Ptr const& digits) + { + + // add the hit to the collection + hits->push_back(hit); + + CreateAssociationsToLastHit(wire, digits); + } // HitCollectionCreator::emplace_back(Hit) + + //---------------------------------------------------------------------- + void HitCollectionCreator::put_into() + { + if (!hits) { + throw art::Exception(art::errors::LogicError) + << "HitCollectionCreator is trying to put into the event" + " a hit collection that was never created!\n"; + } + HitAndAssociationsWriterBase::put_into(); + } // HitCollectionCreator::put_into() + + //---------------------------------------------------------------------- + void HitCollectionCreator::CreateAssociationsToLastHit(art::Ptr const& wire, + art::Ptr const& digits) + { + // if no association is required, we are done + if (!WireAssns && !RawDigitAssns) return; + + // art pointer to the hit we just created + HitPtr_t hit_ptr(CreatePtrToLastHit()); + + // association with wires + if (WireAssns && wire.isNonnull()) + WireAssns->addSingle(wire, hit_ptr); // if it fails, it throws + + // association with wires + if (RawDigitAssns && digits.isNonnull()) + RawDigitAssns->addSingle(digits, hit_ptr); // if it fails, it throws + + } // HitCollectionCreator::CreateAssociationsToLastHit() + + //**************************************************************************** + //*** HitCollectionAssociator + //---------------------------------------------------------------------- + HitCollectionAssociator::HitCollectionAssociator(art::Event& event, + std::string instance_name, + art::InputTag const& WireModuleLabel, + art::InputTag const& RawDigitModuleLabel) + : HitAndAssociationsWriterBase(event, + instance_name, + WireModuleLabel != "", + RawDigitModuleLabel != "") + , wires_label(WireModuleLabel) + , digits_label(RawDigitModuleLabel) + { + hits.reset(new std::vector); + } // HitCollectionAssociator::HitCollectionAssociator() + + //---------------------------------------------------------------------- + sbn::HitCollectionAssociator::HitCollectionAssociator(art::Event& event, + std::string instance_name, + art::InputTag const& WireModuleLabel, + bool doRawDigitAssns /* = false */ + ) + : HitAndAssociationsWriterBase(event, instance_name, WireModuleLabel != "", doRawDigitAssns) + , wires_label(WireModuleLabel) + , digits_label() + { + if (RawDigitAssns && !WireAssns) { + throw art::Exception(art::errors::LogicError) + << "HitCollectionAssociator can't create hit <--> raw digit" + " associations through wires, without wires!\n"; + } + hits.reset(new std::vector); + } // HitCollectionAssociator::HitCollectionAssociator() + + //---------------------------------------------------------------------- + void HitCollectionAssociator::use_hits(std::unique_ptr>&& srchits) + { + hits = std::move(srchits); + } // HitCollectionAssociator::use_hits() + + //---------------------------------------------------------------------- + void HitCollectionAssociator::put_into() + { + prepare_associations(); + HitAndAssociationsWriterBase::put_into(); + } // HitCollectionAssociator::put_into() + + //---------------------------------------------------------------------- + void HitCollectionAssociator::prepare_associations(std::vector const& srchits) + { + if (!RawDigitAssns && !WireAssns) return; // no associations needed + assert(event); + + // we make the associations anew + if (RawDigitAssns) ClearAssociations(*RawDigitAssns); + if (WireAssns) ClearAssociations(*WireAssns); + + // the following is true is we want associations with digits + // but we don't know where digits are; in that case, we try to use wires + const bool bUseWiresForDigits = RawDigitAssns && (digits_label == ""); + + if (WireAssns || bUseWiresForDigits) { + // do we use wires for digit associations too? + + // get the wire collection + art::ValidHandle> hWires = + event->getValidHandle>(wires_label); + + // fill a map of wire index vs. channel number + std::vector WireMap = util::MakeIndex(*hWires, std::mem_fn(&recob::ChannelROI::Channel)); + + // use raw rigit - wire association, assuming they have been produced + // by the same producer as the wire and with the same instance name; + // we don't check whether the data product is found, but the following + // code will have FindOneP throw if that was not the case + // (that's what we would do here anyway, maybe with a better message...) + std::unique_ptr> WireToDigit; + if (bUseWiresForDigits) { + WireToDigit.reset(new art::FindOneP(hWires, *event, wires_label)); + } + + // add associations, hit by hit: + for (size_t iHit = 0; iHit < srchits.size(); ++iHit) { + + // find the channel + size_t iChannel = size_t(srchits[iHit].Channel()); // forcibly converted + + // find the wire associated to that channel + size_t iWire = std::numeric_limits::max(); + if (iChannel < WireMap.size()) iWire = WireMap[iChannel]; + if (iWire == std::numeric_limits::max()) { + throw art::Exception(art::errors::LogicError) + << "No wire associated to channel #" << iChannel << " whence hit #" << iHit + << " comes!\n"; + } // if no channel + + // make the association with wires + if (WireAssns) { + art::Ptr wire(hWires, iWire); + WireAssns->addSingle(wire, CreatePtr(iHit)); + } + + if (bUseWiresForDigits) { + // find the digit associated to that channel + art::Ptr const& digit = WireToDigit->at(iWire); + if (digit.isNull()) { + throw art::Exception(art::errors::LogicError) + << "No raw digit associated to channel #" << iChannel << " whence hit #" << iHit + << " comes!\n"; + } // if no channel + + // make the association + RawDigitAssns->addSingle(digit, CreatePtr(iHit)); + } // if create digit associations through wires + } // for hit + + } // if wire associations + + if (RawDigitAssns && !bUseWiresForDigits) { + // get the digit collection + art::ValidHandle> hDigits = + event->getValidHandle>(digits_label); + + // fill a map of wire index vs. channel number + std::vector DigitMap = + util::MakeIndex(*hDigits, std::mem_fn(&raw::RawDigit::Channel)); + + // add associations, hit by hit: + for (size_t iHit = 0; iHit < srchits.size(); ++iHit) { + + // find the channel + size_t iChannel = size_t(srchits[iHit].Channel()); // forcibly converted + + // find the digit associated to that channel + size_t iDigit = std::numeric_limits::max(); + if (iChannel < DigitMap.size()) iDigit = DigitMap[iChannel]; + if (iDigit == std::numeric_limits::max()) { + throw art::Exception(art::errors::LogicError) + << "No raw digit associated to channel #" << iChannel << " whence hit #" << iHit + << " comes!\n"; + } // if no channel + + // make the association + art::Ptr digit(hDigits, iDigit); + RawDigitAssns->addSingle(digit, CreatePtr(iHit)); + + } // for hit + } // if we have rawdigit label + + } // HitCollectionAssociator::put_into() + + //**************************************************************************** + //*** HitRefinerAssociator + //---------------------------------------------------------------------- + HitRefinerAssociator::HitRefinerAssociator(art::Event& event, + art::InputTag const& HitModuleLabel, + std::string instance_name /* = "" */, + bool doWireAssns /* = true */, + bool doRawDigitAssns /* = true */ + ) + : HitAndAssociationsWriterBase(event, instance_name, doWireAssns, doRawDigitAssns) + , hits_label(HitModuleLabel) + { + hits.reset(new std::vector); + } // HitRefinerAssociator::HitRefinerAssociator() + + //---------------------------------------------------------------------- + void HitRefinerAssociator::use_hits(std::unique_ptr>&& srchits) + { + hits = std::move(srchits); + } // HitRefinerAssociator::use_hits() + + //---------------------------------------------------------------------- + void HitRefinerAssociator::put_into() + { + prepare_associations(); + HitAndAssociationsWriterBase::put_into(); + } // HitRefinerAssociator::put_into() + + //---------------------------------------------------------------------- + void HitRefinerAssociator::prepare_associations(std::vector const& srchits) + { + if (!RawDigitAssns && !WireAssns) return; // no associations needed + assert(event); + + // we make the associations anew + if (RawDigitAssns) ClearAssociations(*RawDigitAssns); + + // read the hits; this is going to hurt performances... + // no solution to that until there is a way to have a lazy read + art::ValidHandle> hHits = + event->getValidHandle>(hits_label); + + // now get the associations + if (WireAssns) { + // we make the associations anew + ClearAssociations(*WireAssns); + + // find the associations between the hits and the wires + art::FindOneP HitToWire(hHits, *event, hits_label); + if (!HitToWire.isValid()) { + throw art::Exception(art::errors::ProductNotFound) + << "Can't find the associations between hits and wires produced by '" << hits_label + << "'!\n"; + } // if no association + + // fill a map of wire vs. channel number + std::vector> WireMap; + for (size_t iAssn = 0; iAssn < HitToWire.size(); ++iAssn) { + art::Ptr wire = HitToWire.at(iAssn); + if (wire.isNull()) continue; + size_t channelID = (size_t)wire->Channel(); + if (WireMap.size() <= channelID) // expand the map of necessary + WireMap.resize(std::max(channelID + 1, 2 * WireMap.size()), {}); + WireMap[channelID] = std::move(wire); + } // for + + // now go through all the hits... + for (size_t iHit = 0; iHit < srchits.size(); ++iHit) { + recob::Hit const& hit = srchits[iHit]; + size_t channelID = (size_t)hit.Channel(); + + // no association if there is no wire to associate with + if ((channelID >= WireMap.size()) || !WireMap[channelID]) continue; + + // create an association using the same wire pointer + WireAssns->addSingle(WireMap[channelID], CreatePtr(iHit)); + } // for hits + } // if wire associations + + // now get the associations + if (RawDigitAssns) { + // we make the associations anew + ClearAssociations(*RawDigitAssns); + + // find the associations between the hits and the raw digits + art::FindOneP HitToDigits(hHits, *event, hits_label); + if (!HitToDigits.isValid()) { + throw art::Exception(art::errors::ProductNotFound) + << "Can't find the associations between hits and raw digits" + << " produced by '" << hits_label << "'!\n"; + } // if no association + + // fill a map of digits vs. channel number + std::vector> DigitMap; + for (size_t iAssn = 0; iAssn < HitToDigits.size(); ++iAssn) { + art::Ptr digits = HitToDigits.at(iAssn); + if (digits.isNull()) continue; + size_t channelID = (size_t)digits->Channel(); + if (DigitMap.size() <= channelID) // expand the map of necessary + DigitMap.resize(std::max(channelID + 1, 2 * DigitMap.size()), {}); + DigitMap[channelID] = std::move(digits); + } // for + + // now go through all the hits... + for (size_t iHit = 0; iHit < srchits.size(); ++iHit) { + recob::Hit const& hit = srchits[iHit]; + size_t channelID = (size_t)hit.Channel(); + + // no association if there is no digits to associate with + if ((channelID >= DigitMap.size()) || !DigitMap[channelID]) continue; + + // create an association using the same digits pointer + RawDigitAssns->addSingle(DigitMap[channelID], CreatePtr(iHit)); + } // for hits + } // if digit associations + + } // HitRefinerAssociator::put_into() + +} // namespace recob diff --git a/sbncode/HitFinder/HitFinderUtilities/HitCreator.h b/sbncode/HitFinder/HitFinderUtilities/HitCreator.h new file mode 100644 index 000000000..658fc1a9b --- /dev/null +++ b/sbncode/HitFinder/HitFinderUtilities/HitCreator.h @@ -0,0 +1,1062 @@ +/** **************************************************************************** + * @file HitCreator.h + * @brief Helper functions to create a hit from ChannelROIs + * @date April 20, 2025 + * @author Tracy Usher (usher@slac.stanford.edu) + * @see Hit.h HitCreator.cxx + * + * ****************************************************************************/ + +#ifndef SBN_ARTDATAHELPERS_HITCREATOR_H +#define SBN_ARTDATAHELPERS_HITCREATOR_H + +// LArSoft libraries +#include "lardataobj/RawData/RawDigit.h" +#include "lardataobj/RecoBase/Hit.h" +#include "sbnobj/ICARUS/TPC/ChannelROI.h" + +// framework libraries +#include "art/Framework/Core/ProducesCollector.h" +#include "art/Persistency/Common/PtrMaker.h" +#include "canvas/Persistency/Common/Assns.h" +#include "canvas/Persistency/Common/Ptr.h" +#include "canvas/Utilities/Exception.h" +#include "canvas/Utilities/InputTag.h" + +// C/C++ standard library +#include +#include // std::move() +#include + +namespace geo { + struct WireID; +} +namespace raw { + class RawDigit; +} +namespace art { + class ProducesCollector; + class Event; +} + +/// Reconstruction base classes +namespace sbn { + + /** ************************************************************************** + * @brief Class managing the creation of a new `recob::Hit` object. + * + * In order to be as simple as possible (Plain Old Data), data products like + * `recob::Hit` need to be stripped of most of their functions, including the + * ability to communicate whether a value we try to store is invalid + * (that would require a art::Exception` -- art -- or at least a message on + * the screen -- MessageFacility) and the ability to read things from event, + * services (e.g. geometry) etc. + * + * A Creator is a class that creates a temporary data product, and at the + * end it yields it to the caller for storage. + * This last step should be by move construction, although a copy method is + * also provided. + * + * An example of creating a `recob::Hit` object (assuming all the relevant + * variables have been assigned proper values): + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * recob::HitCreator hit( + * wire, wireID, + * start_tick, end_tick, rms, + * peak_time, sigma_peak_time, peak_amplitude, sigma_peak_amplitude, + * hit_integral, hit_sigma_integral, summedADC, + * multiplicity, local_index, goodness_of_fit, dof + * ); + * hit.push_back(hit.move()); // hit content is not valid any more + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * + * This is a one-step creation object: the hit is constructed at the same + * time the HitCreator is, and no facility is offered to modify the + * constructed hit, or to create another one. + * + * The constructors currently provided are: + * 1. from RawDigit (extracts channel, view and signal type [CVS] thanks to + * geometry) + * 2. from `recob::ChannelROI`, [CVS] + * 3. from `recob::ChannelROI`, [CVS], `summedADC` is automatically computed from + * wire + * 4. from `recob::ChannelROI`, [CVS], start and stop time from a region of interest + * 5. from `recob::ChannelROI`, [CVS], start and stop time from index of region of + * interest + */ + class HitCreator { + public: + /// Type of one region of interest. + using RegionOfInterest_t = recob::ChannelROI::RegionsOfInterest_f::datarange_t; + + // destructor, copy and move constructor and assignment as default + + /** + * @brief Constructor: extracts some information from raw digit. + * @param digits a pointer to a `raw::RawDigit` (for channel, view, signal + * type) + * @param wireID ID of the wire the hit is on + * @param start_tick first tick in the region the hit was extracted from + * @param end_tick first tick after the region the hit was extracted from + * @param rms RMS of the signal hit, in TDC time units + * @param peak_time time at peak of the signal, in TDC time units + * @param sigma_peak_time uncertainty on time at peak, in TDC time units + * @param peak_amplitude amplitude of the signal at peak, in ADC units + * @param sigma_peak_amplitude uncertainty on amplitude at peak + * @param hit_integral total charge integrated under the hit signal + * @param hit_sigma_integral uncertainty on the total hit charge + * @param ROIsummedADC total ADC count in the region assigned to the hit ROI + * @param HitsummedADC total ADC count in the region assigned to the hit HIT + * @param multiplicity number of hits in the region it was extracted from + * @param local_index index of this hit in the region it was extracted + * from + * @param goodness_of_fit quality parameter for the hit + * @param dof degrees of freedom in the definition of the hit shape + * + * The information used from the raw digit is the channel ID; view and + * signal type are obtained from geometry. + */ + HitCreator(raw::RawDigit const& digits, + geo::WireID const& wireID, + raw::TDCtick_t start_tick, + raw::TDCtick_t end_tick, + float rms, + float peak_time, + float sigma_peak_time, + float peak_amplitude, + float sigma_peak_amplitude, + float hit_integral, + float hit_sigma_integral, + float ROIsummedADC, + float HitsummedADC, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof); + + /** + * @brief Constructor: extracts some information from wire. + * @param wire a pointer to a `recob::ChannelROI` (for channel, view, signal + * type) + * @param wireID ID of the wire the hit is on + * @param start_tick first tick in the region the hit was extracted from + * @param end_tick first tick after the region the hit was extracted from + * @param rms RMS of the signal hit, in TDC time units + * @param peak_time time at peak of the signal, in TDC time units + * @param sigma_peak_time uncertainty on time at peak, in TDC time units + * @param peak_amplitude amplitude of the signal at peak, in ADC units + * @param sigma_peak_amplitude uncertainty on amplitude at peak + * @param hit_integral total charge integrated under the hit signal + * @param hit_sigma_integral uncertainty on the total hit charge + * @param ROIsummedADC total ADC count in the region assigned to the hit ROI + * @param HitsummedADC total ADC count in the region assigned to the hit Hit + * @param multiplicity number of hits in the region it was extracted from + * @param local_index index of this hit in the region it was extracted + * from + * @param goodness_of_fit quality parameter for the hit + * @param dof degrees of freedom in the definition of the hit shape + * + * The information used from the wire are the channel ID and view; + * the signal type is obtained from geometry. + */ + HitCreator(recob::ChannelROI const& wire, + geo::WireID const& wireID, + raw::TDCtick_t start_tick, + raw::TDCtick_t end_tick, + float rms, + float peak_time, + float sigma_peak_time, + float peak_amplitude, + float sigma_peak_amplitude, + float hit_integral, + float hit_sigma_integral, + float ROIsummedADC, + float HitsummedADC, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof); + + /** + * @brief Constructor: computes sum of ADC from wire. + * @param wire a pointer to a `recob::ChannelROI` (for channel, view, signal + * type) + * @param wireID ID of the wire the hit is on + * @param start_tick first tick in the region the hit was extracted from + * @param end_tick first tick after the region the hit was extracted from + * @param rms RMS of the signal hit, in TDC time units + * @param peak_time time at peak of the signal, in TDC time units + * @param sigma_peak_time uncertainty on time at peak, in TDC time units + * @param peak_amplitude amplitude of the signal at peak, in ADC units + * @param sigma_peak_amplitude uncertainty on amplitude at peak + * @param hit_integral total charge integrated under the hit signal + * @param hit_sigma_integral uncertainty on the total hit charge + * @param multiplicity number of hits in the region it was extracted from + * @param local_index index of this hit in the region it was extracted from + * @param goodness_of_fit quality parameter for the hit + * @param dof degrees of freedom in the definition of the hit shape + * + * The information used from the wire are the channel ID, view; + * the signal type is obtained from geometry. + * + * The sum of ADC counts is automatically computed over the whole range + * of the wire signal between `start_tick` and `end_tick` + * (the latter excluded). + */ + HitCreator(recob::ChannelROI const& wire, + geo::WireID const& wireID, + raw::TDCtick_t start_tick, + raw::TDCtick_t end_tick, + float rms, + float peak_time, + float sigma_peak_time, + float peak_amplitude, + float sigma_peak_amplitude, + float hit_integral, + float hit_sigma_integral, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof); + + /** + * @brief Constructor: uses region of interest specified by index. + * @param wire a pointer to a `recob::ChannelROI` (for channel, view, signal + * type) + * @param wireID ID of the wire the hit is on + * @param rms RMS of the signal hit, in TDC time units + * @param peak_time time at peak of the signal, in TDC time units + * @param sigma_peak_time uncertainty on time at peak, in TDC time units + * @param peak_amplitude amplitude of the signal at peak, in ADC units + * @param sigma_peak_amplitude uncertainty on amplitude at peak + * @param hit_integral total charge integrated under the hit signal + * @param hit_sigma_integral uncertainty on the total hit charge + * @param ROIsummedADC total ADC count in the region assigned to the hit ROI + * @param HitsummedADC total ADC count in the region assigned to the hit Hit + * @param multiplicity number of hits in the region it was extracted from + * @param local_index index of this hit in the region it was extracted + * from + * @param goodness_of_fit quality parameter for the hit + * @param dof degrees of freedom in the definition of the hit shape + * @param signal the signal region the hit was extracted from + * + * The information used from the wire are the channel ID, view + * and the region of interest; the signal type is obtained from + * geometry. + * + * Signal start and end ticks are extracted from the region of interest. + */ + HitCreator(recob::ChannelROI const& wire, + geo::WireID const& wireID, + float rms, + float peak_time, + float sigma_peak_time, + float peak_amplitude, + float sigma_peak_amplitude, + float hit_integral, + float hit_sigma_integral, + float ROIsummedADC, + float HitsummedADC, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof, + RegionOfInterest_t const& signal); + + /** + * @brief Constructor: uses region of interest specified by index. + * @param wire a pointer to a `recob::ChannelROI` (for channel, view, signal + * type) + * @param wireID ID of the wire the hit is on + * @param rms RMS of the signal hit, in TDC time units + * @param peak_time time at peak of the signal, in TDC time units + * @param sigma_peak_time uncertainty on time at peak, in TDC time units + * @param peak_amplitude amplitude of the signal at peak, in ADC units + * @param sigma_peak_amplitude uncertainty on amplitude at peak + * @param hit_integral total charge integrated under the hit signal + * @param hit_sigma_integral uncertainty on the total hit charge + * @param ROI summedADC total ADC count in the region assigned to the hit ROI + * @param Hit summedADC total ADC count in the region assigned to the hit Hit + * @param multiplicity number of hits in the region it was extracted from + * @param local_index index of this hit in the region it was extracted from + * @param goodness_of_fit quality parameter for the hit + * @param dof degrees of freedom in the definition of the hit shape + * @param iSignalRoI index in the wire of the signal region the hit was + * extracted from + * + * The information used from the wire are the channel ID, view + * and the region of interest; the signal type is obtained from + * geometry. + * + * Signal start and end ticks are extracted from the region of interest. + */ + HitCreator(recob::ChannelROI const& wire, + geo::WireID const& wireID, + float rms, + float peak_time, + float sigma_peak_time, + float peak_amplitude, + float sigma_peak_amplitude, + float hit_integral, + float hit_sigma_integral, + float ROIsummedADC, + float HitsummedADC, + short int multiplicity, + short int local_index, + float goodness_of_fit, + int dof, + size_t iSignalRoI); + + /** + * @brief Constructor: copies from an existing hit. + * @param from the original hit + */ + HitCreator(recob::Hit const& from); + + /** + * @brief Constructor: copies from an existing hit, changing wire ID. + * @param from the original hit + * @param wireID ID of the new wire the hit is on + */ + HitCreator(recob::Hit const& from, geo::WireID const& wireID); + + /** + * @brief Prepares the constructed hit to be moved away. + * @return a right-value reference to the constructed hit + * + * Despite the name, no move happens in this function. + * Move takes place in the caller code as proper; for example: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * // be hit a HitCreator instance: + * std::vector Hits; + * hit.move(); // nothing happens + * Hits.push_back(hit.move()); // here the copy happens + * recob::Hit single_hit(hit.move()); // wrong! hit is empty now + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * + */ + recob::Hit&& move() { return std::move(hit); } + + /** + * @brief Returns the constructed wire + * @return a constant reference to the constructed wire + * + * Despite the name, no copy happens in this function. + * Copy takes place in the caller code as proper; for example: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * // be Hit a HitCreator instance: + * std::vector Hits; + * hit.copy(); // nothing happens + * Hits.push_back(hit.copy()); // here a copy happens + * recob::Hit single_hit(hit.copy()); // hit is copied again + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * + */ + recob::Hit const& copy() const { return hit; } + + protected: + recob::Hit hit; ///< Local instance of the hit being constructed. + + }; // class HitCreator + + /** ************************************************************************** + * @brief Base class handling a collection of hits and its associations. + * + * Instead of creating a collection of hits, one for its association with + * wires and one for its association with raw digits, one can use a class + * derived from this one: + * - `HitCollectionCreator`: push new hits one by one + * - `HitCollectionAssociator`: push a complete collection of hits + * - `HitRefinerAssociator`: push a complete collection of hits deriving their + * associations from other hits + * Using `put_into()` will transfer into the event the data. + * + * The typical usage is to have the constructor of the module call the static + * function + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * HitAndAssociationsWriterBase::declare_products(*this); + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * + * (this example declares a collection with empty instance name and that we + * want associations to both wires and raw digits), and then in `produce()`: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * HitAndAssociationsWriterDerived hcol(*this, event); + * + * // ... fill hcol in the proper way ... + * + * hcol.put_into(); // calls art::Event::put() + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * + */ + class HitAndAssociationsWriterBase { + public: + // no public constructor: use one of the derived classes! + // destructor, copy and move constructors and assignment are default + + /// Returns the number of hits currently in the collection. + size_t size() const { return hits ? hits->size() : 0; } + + /** + * @brief Moves the data into the event. + * + * The calling module must have already declared the production of these + * products with the proper instance name. + * After the move, the collections in this object are empty. + * + * @deprecated Use the version with no arguments instead. + */ + void put_into(art::Event&) { put_into(); } + + /** + * @brief Moves the data into the event. + * + * The calling module must have already declared the production of these + * products with the proper instance name. + * After the move, the collections in this object are empty. + */ + void put_into(); + + /// Returns a read-only reference to the current list of hits. + std::vector const& peek() const { return *hits; } + + /** + * @brief Declares the hit products we are going to fill. + * @tparam ModuleType type of producing module (`EDProducer` or `EDFilter`) + * @param producer the module producing the data products + * @param instance_name name of the instance for all data products + * @param doWireAssns whether to enable associations to wires + * @param doRawDigitAssns whether to enable associations to raw digits + * + * This declaration must be given in the constructor of producer. + * It is equivalent to manually declare the relevant among these products: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * produces>(prod_instance); + * produces>(prod_instance); + * produces>(prod_instance); + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * in the producer constructor. + * All the data products (hit collection and associations) will have the + * specified product instance name. + */ + static void declare_products(art::ProducesCollector& collector, + std::string instance_name = "", + bool doWireAssns = true, + bool doRawDigitAssns = true); + + protected: + using HitPtr_t = art::Ptr; ///< Type of art pointer to Hit. + + std::string prod_instance; ///< Tame of the instance for data products. + + /// Collection of hits. + std::unique_ptr> hits; + /// Associations with wires. + std::unique_ptr> WireAssns; + /// Associations with raw digits. + std::unique_ptr> RawDigitAssns; + + art::Event* event = nullptr; ///< Pointer to the event we are using. + + art::PtrMaker hitPtrMaker; ///< Tool to create hit pointers, + + /** + * @brief Constructor: sets instance name and whether to build associations. + * @param event the event the products are going to be put into + * @param instance_name name of the instance for all data products + * @param doWireAssns whether to enable associations to wires + * @param doRawDigitAssns whether to enable associations to raw digits + * + * All the data products (hit collection and associations) will have the + * specified product instance name. + */ + HitAndAssociationsWriterBase(art::Event& event, + std::string instance_name, + bool doWireAssns, + bool doRawDigitAssns); + + /// Creates an art pointer to the hit with the specified index. + HitPtr_t CreatePtr(size_t index) const { return hitPtrMaker(index); } + + }; // class HitAndAssociationsWriterBase + + /** ************************************************************************** + * @brief A class handling a collection of hits and its associations. + * + * Instead of creating a collection of hits, one for its association with + * wires and one for its association with raw digits, one can push hits into + * this object, and then move it into the event. + */ + class HitCollectionCreator : public HitAndAssociationsWriterBase { + public: + /// @name Constructors + /// @{ + /** + * @brief Constructor: sets instance name and whether to build associations. + * @param event the event the products are going to be put into + * @param instance_name name of the instance for all data products + * @param doWireAssns whether to enable associations to wires + * @param doRawDigitAssns whether to enable associations to raw digits + * + * All the data products (hit collection and associations) will have the + * specified product instance name. + */ + HitCollectionCreator(art::Event& event, + std::string instance_name = "", + bool doWireAssns = true, + bool doRawDigitAssns = true); + + /** + * @brief Constructor: no product instance name. + * @param event the event the products are going to be put into + * @param doWireAssns whether to enable associations to wires + * @param doRawDigitAssns whether to enable associations to raw digits + */ + HitCollectionCreator(art::Event& event, bool doWireAssns, bool doRawDigitAssns) + : HitCollectionCreator(event, "", doWireAssns, doRawDigitAssns) + {} + + /// @} + + // destructor, copy and move constructors and assignment are default + + /// @name Addition of hits + /// @{ + /** + * @brief Adds the specified hit to the data collection. + * @param hit the hit that will be moved into the collection + * @param wire art pointer to the wire to be associated to this hit + * @param digits art pointer to the raw digits to be associated to this hit + * + * After this call, hit will be invalid. + * If a art pointer is not valid, that association will not be stored. + */ + void emplace_back(recob::Hit&& hit, + art::Ptr const& wire = art::Ptr(), + art::Ptr const& digits = art::Ptr()); + + /** + * @brief Adds the specified hit to the data collection. + * @param hit the hit that will be copied into the collection + * @param wire art pointer to the wire to be associated to this hit + * @param digits art pointer to the raw digits to be associated to this hit + * + * If a art pointer is not valid, that association will not be stored. + */ + void emplace_back(recob::Hit const& hit, + art::Ptr const& wire = art::Ptr(), + art::Ptr const& digits = art::Ptr()); + + /** + * @brief Adds the specified hit to the data collection. + * @param hit the HitCreator object containing the hit + * @param wire art pointer to the wire to be associated to this hit + * @param digits art pointer to the raw digits to be associated to this hit + * + * After this call, the hit creator will be empty. + * If a art pointer is not valid, that association will not be stored. + */ + void emplace_back(HitCreator&& hit, + art::Ptr const& wire = art::Ptr(), + art::Ptr const& digits = art::Ptr()) + { + emplace_back(hit.move(), wire, digits); + } + + /** + * @brief Adds the specified hit to the data collection. + * @param hit the hit that will be moved into the collection + * @param digits art pointer to the raw digits to be associated to this hit + * + * After this call, hit will be invalid. + * If the digit pointer is not valid, its association will not be stored. + */ + void emplace_back(recob::Hit&& hit, art::Ptr const& digits) + { + emplace_back(std::move(hit), art::Ptr(), digits); + } + + /** + * @brief Adds the specified hit to the data collection. + * @param hit the HitCreator object containing the hit + * @param digits art pointer to the raw digits to be associated to this hit + * + * After this call, the hit creator will be empty. + * If the digit pointer is not valid, its association will not be stored. + */ + void emplace_back(HitCreator&& hit, art::Ptr const& digits) + { + emplace_back(std::move(hit), art::Ptr(), digits); + } + + /** + * @brief Adds the specified hit to the data collection. + * @param hit the HitCreator object containing the hit + * @param digits art pointer to the raw digits to be associated to this hit + * + * If the digit pointer is not valid, its association will not be stored. + */ + void emplace_back(HitCreator const& hit, art::Ptr const& digits) + { + emplace_back(std::move(hit.copy()), art::Ptr(), digits); + } + /// @} + + /// Returns the number of hits currently in the collection. + size_t size() const { return hits->size(); } + + /// Prepares the collection to host at least `new_size` hits. + void reserve(size_t new_size) + { + if (hits) hits->reserve(new_size); + } + + /** + * @brief Moves the data into an event. + * + * The calling module must have already declared the production of these + * products with the proper instance name. + * After the move, the collections in this object are empty. + * + * @deprecated Use the version with no arguments instead. + */ + void put_into(art::Event&) { put_into(); } + + /** + * @brief Moves the data into the event. + * + * The calling module must have already declared the production of these + * products with the proper instance name. + * After the move, the collections in this object are empty. + */ + void put_into(); + + /// Returns a read-only reference to the current list of hits. + std::vector const& peek() const { return *hits; } + + protected: + using HitPtr_t = HitAndAssociationsWriterBase::HitPtr_t; + + /// Creates an art pointer to the hit with the last index. + HitPtr_t CreatePtrToLastHit() const + { + return hits->empty() ? HitPtr_t() : CreatePtr(hits->size() - 1); + } + + /// Creates associations between the last hit and the specified pointers. + void CreateAssociationsToLastHit(art::Ptr const& wire, + art::Ptr const& digits); + + }; // class HitCollectionCreator + + /** ************************************************************************** + * @brief A class handling a collection of hits and its associations. + * + * Use this object if you already have a collection of `recob::Hit` and you + * simply want the hits associated to the wire and digit with the same + * channel. + */ + class HitCollectionAssociator : public HitAndAssociationsWriterBase { + public: + /// @name Constructors + /// @{ + /** + * @brief Constructor: sets instance name and whether to build associations. + * @param event the event the products are going to be put into + * @param instance_name name of the instance for all data products + * @param WireModuleLabel label of the module used to create wires + * @param RawDigitModuleLabel label of the module used to create raw digits + * + * All the data products (hit collection and associations) will have the + * specified product instance name. + * + * If a label is empty, the corresponding association will not be produced. + */ + HitCollectionAssociator(art::Event& event, + std::string instance_name, + art::InputTag const& WireModuleLabel, + art::InputTag const& RawDigitModuleLabel); + + /** + * @brief Constructor: sets instance name and whether to build associations. + * @param event the event the products are going to be put into + * @param WireModuleLabel label of the module used to create wires + * @param RawDigitModuleLabel label of the module used to create raw digits + * + * All the data products (hit collection and associations) will have a + * default, empty product instance name. + * + * If a label is empty, the corresponding association will not be produced. + */ + HitCollectionAssociator(art::Event& event, + art::InputTag const& WireModuleLabel, + art::InputTag const& RawDigitModuleLabel) + : HitCollectionAssociator(event, "", WireModuleLabel, RawDigitModuleLabel) + {} + + /** + * @brief Constructor: sets instance name and whether to build associations. + * @param event the event the products are going to be put into + * @param instance_name name of the instance for all data products + * @param WireModuleLabel label of the module used to create wires + * @param doRawDigitAssns whether to write associations with raw digits + * + * All the data products (hit collection and associations) will have the + * specified product instance name. + * + * The raw digit association is built out of their existing associations + * with wires, rather than by directly using the raw digits data product. + */ + HitCollectionAssociator(art::Event& event, + std::string instance_name, + art::InputTag const& WireModuleLabel, + bool doRawDigitAssns); + + /** + * @brief Constructor: sets instance name and whether to build associations. + * @param event the event the products are going to be put into + * @param WireModuleLabel label of the module used to create wires + * @param doRawDigitAssns whether to write associations with raw digits + * + * All the data products (hit collection and associations) will have the + * default, empty product instance name. + * + * The raw digit association is built out of their existing associations + * with wires, rather than by directly using the raw digits data product. + */ + HitCollectionAssociator(art::Event& event, + art::InputTag const& WireModuleLabel, + bool doRawDigitAssns) + : HitCollectionAssociator(event, "", WireModuleLabel, doRawDigitAssns) + {} + + /// @} + + // destructor, copy and move constructors and assignment are default + + /** + * @brief Uses the specified collection as data product. + * @param srchits the collection to be used as data product + * + * The very same collection is put into the event. + * This object will temporary own the collection until the hits are put into + * the event. + * If there were previous hits in the object, they are lost. + */ + void use_hits(std::unique_ptr>&& srchits); + + /** + * @brief Moves the data into the event. + * + * The calling module must have already declared the production of these + * products with the proper instance name. + * After the move, the collections in this object are empty. + * + * @deprecated Use the version with no arguments instead. + */ + void put_into(art::Event&) { put_into(); } + + /** + * @brief Moves the data into the event. + * + * The calling module must have already declared the production of these + * products with the proper instance name. + * After the move, the collections in this object are empty. + */ + void put_into(); + + protected: + /// Label of the collection of wires to associate. + art::InputTag wires_label; + /// Label of raw digits collection to associate. + art::InputTag digits_label; + + /// Finds out the associations for the specified hits. + void prepare_associations(std::vector const& srchits); + + /// Finds out the associations for the current hits. + void prepare_associations() { prepare_associations(*hits); } + + }; // class HitCollectionAssociator + + /** ************************************************************************** + * @brief A class handling a collection of hits and its associations. + * + * Use this object if you already have a `recob::Hit` data product and + * another collection that is going to become a data product, and you + * simply want the new hits associated to the wire and digit with the same + * channel. + * No hit-to-hit association is attempted (that would be, incidentally, not + * supported by art): the data product is used to get all the associated + * wires and digits, then they are associated to the new hits by channel ID. + * If a channel is not available, a warning is produced. If different hits + * on the same channel are associated to different wires or raw digits, an + * exception is thrown. + */ + class HitRefinerAssociator : public HitAndAssociationsWriterBase { + public: + /// @name Constructors + /// @{ + /** + * @brief Constructor: sets instance name and whether to build associations. + * @param event the event the products are going to be put into + * @param HitModuleLabel label of the module used to create hits + * @param instance_name name of the instance for all data products + * @param doWireAssns whether to enable associations to wires + * @param doRawDigitAssns whether to enable associations to raw digits + * + * All the data products (hit collection and associations) will have the + * specified product instance name. + */ + HitRefinerAssociator(art::Event& event, + art::InputTag const& HitModuleLabel, + std::string instance_name = "", + bool doWireAssns = true, + bool doRawDigitAssns = true); + + /** + * @brief Constructor: sets instance name and whether to build associations. + * @param event the event the products are going to be put into + * @param HitModuleLabel label of the module used to create hits + * @param doWireAssns whether to enable associations to wires + * @param doRawDigitAssns whether to enable associations to raw digits + * + * All the data products (hit collection and associations) will have an + * empty product instance name. + */ + HitRefinerAssociator(art::Event& event, + art::InputTag const& HitModuleLabel, + bool doWireAssns, + bool doRawDigitAssns = true) + : HitRefinerAssociator(event, HitModuleLabel, "", doWireAssns, doRawDigitAssns) + {} + + /// @} + + // destructor, copy and move constructors and assignment are default + + /** + * @brief Uses the specified collection as data product. + * @param srchits the collection to be used as data product + * + * The very same collection is put into the event. + * This object will temporary own the collection until the hits are put into + * the event. + * If there were previous hits in the object, they are lost. + */ + void use_hits(std::unique_ptr>&& srchits); + + /** + * @brief Moves the data into the event. + * + * The calling module must have already declared the production of these + * products with the proper instance name. + * After the move, the collections in this object are empty. + * + * @deprecated Use the version with no arguments instead. + * + */ + void put_into(art::Event&) { put_into(); } + + /** + * @brief Moves the data into the event. + * + * The calling module must have already declared the production of these + * products with the proper instance name. + * After the move, the collections in this object are empty. + */ + void put_into(); + + protected: + art::InputTag hits_label; ///< Label of the collection of hits. + + /// Finds out the associations for the specified hits. + void prepare_associations(std::vector const& srchits); + + /// Finds out the associations for the current hits. + void prepare_associations() { prepare_associations(*hits); } + + }; // class HitRefinerAssociator + + // --------------------------------------------------------------------------- + /** + * @brief A helper to centralise creation of a hit collection data product. + * @tparam Writer writer class to manage + * @tparam ModuleType owning module: `art::EDProducer` or `art::EDFilter` + * + * This class adds an indirection layer to the model relying on + * `HitAndAssociationsWriter`. In that one, two different steps are required, + * one in the constructor of the module, where data products are declared, and + * one in the place where hits are actually assembled. + * These two steps need consistent setup, but they are separate and + * formally independent. The "manager" approach consists of an object + * performing the first step directly, and delivering an already configured + * object for the second step. + * + * An example of usage in a module: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * class MyHitProducer: public art::EDProducer { + * + * recob::HitAndAssociationsWriterManager + * hitCollCreator; + * + * public: + * + * MyHitProducer(fhicl::ParameterSet const& pset) + * : EDProducer{pset} + * , hitCollCreator(*this, pset.get("instanceName", "")) + * {} + * + * void produce(art::Event& event) + * { + * auto hitCollWriter = hitCollCreator.collectionWriter(event); + * + * for (recob::ChannelROI const& wire: Wires) { + * // create hits... + * hitCollWriter.emplace_back(hit, wire, digit); + * } + * hitCollWriter.put_into(); + * } + * + * }; // class MyHitProducer + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * + */ + template + class HitAndAssociationsWriterManager { + + public: + using Writer_t = Writer; ///< Type of managed hit collection writer. + + /** + * @brief Constructor: does not declare anything. + * + * This constructor does not declare products. Calling `declare_products()` + * explicitly is then required in the module constructor. + * + */ + HitAndAssociationsWriterManager() = default; + + /** + * @brief Declares the hit products we are going to fill. + * @param collector the module this manager is bound to + * @param instanceName name of the instance for all data products + * @param doWireAssns whether to enable associations to wires + * @param doRawDigitAssns whether to enable associations to raw digits + * + * This constructor calls `declareProducts()`. + */ + HitAndAssociationsWriterManager(art::ProducesCollector& collector, + std::string instanceName = "", + bool doWireAssns = true, + bool doRawDigitAssns = true); + + /** + * @brief Declares the hit products we are going to fill. + * @param collector the module this manager is bound to + * @param instanceName name of the instance for all data products + * @param doWireAssns whether to enable associations to wires + * @param doRawDigitAssns whether to enable associations to raw digits + * + * This declaration must be made in the constructor of producer. + * It is equivalent to manually declare the relevant among these products: + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{.cpp} + * produces>(prod_instance); + * produces>(prod_instance); + * produces>(prod_instance); + * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + * in the producer constructor. + * All the data products (hit collection and associations) will have the + * specified product instance name. + */ + void declareProducts(art::ProducesCollector& collector, + std::string instanceName = "", + bool doWireAssns = true, + bool doRawDigitAssns = true); + + /// Returns a new writer already configured. + Writer_t collectionWriter(art::Event& event) const; + + /// Returns the configured product instance name. + std::string instanceName() const { return prodInstance; } + + /// Returns whether the class is fully configured. + bool ready() const noexcept { return collector_p != nullptr; } + + protected: + art::ProducesCollector* collector_p = nullptr; ///< Collector this manager is bound to. + + std::string prodInstance; ///< Tame of the instance for data products. + + /// Whether we produce hit-digit associations. + bool hasRawDigitAssns = true; + + /// Whether we produce hit-wire associations. + bool hasWireAssns = true; + + }; // class HitAndAssociationsWriterManager + + /// A manager for `recob::HitCollectionCreator` writer class. + using HitCollectionCreatorManager = HitAndAssociationsWriterManager; + +} // namespace sbn + +//------------------------------------------------------------------------------ +//--- template implementation +//------------------------------------------------------------------------------ +//--- recob::HitAndAssociationsWriterBase +//--- +//------------------------------------------------------------------------------ +//--- recob::HitAndAssociationsWriterManager +//--- +template +sbn::HitAndAssociationsWriterManager::HitAndAssociationsWriterManager( + art::ProducesCollector& collector, + std::string instanceName /* = "" */, + bool doWireAssns /* = true */, + bool doRawDigitAssns /* = true */ +) +{ + declareProducts(collector, instanceName, doWireAssns, doRawDigitAssns); +} // recob::HitAndAssociationsWriterManager::HitAndAssociationsWriterManager() + +//------------------------------------------------------------------------------ +template +void sbn::HitAndAssociationsWriterManager::declareProducts( + art::ProducesCollector& collector, + std::string instanceName /* = "" */, + bool doWireAssns /* = true */, + bool doRawDigitAssns /* = true */ +) +{ + if (collector_p) { + // this means you already called to declaredProducts() + // or used the wrong constructor (which did that for you): + throw art::Exception(art::errors::LogicError) + << "HitAndAssociationsWriter<> has already declared its products."; + } + collector_p = &collector; + prodInstance = instanceName; + hasWireAssns = doWireAssns; + hasRawDigitAssns = doRawDigitAssns; + HitAndAssociationsWriterBase::declare_products( + collector, prodInstance, hasWireAssns, hasRawDigitAssns); +} // recob::HitAndAssociationsWriterManager::declareProducts() + +//------------------------------------------------------------------------------ +template +typename sbn::HitAndAssociationsWriterManager::Writer_t +sbn::HitAndAssociationsWriterManager::collectionWriter(art::Event& event) const +{ + if (!collector_p) { + // this means you forgot to code a call to declaredProducts() + // or used the wrong constructor: + throw art::Exception(art::errors::LogicError) + << "HitAndAssociationsWriter<>::collectionWriter() called" + " before products are declared."; + } + return {event, prodInstance, hasWireAssns, hasRawDigitAssns}; +} // recob::HitAndAssociationsWriterManager::collectionWriter() + +//------------------------------------------------------------------------------ + +#endif // LARDATA_ARTDATAHELPERS_HITCREATOR_H diff --git a/sbncode/HitFinder/WireToChannelROI_module.cc b/sbncode/HitFinder/WireToChannelROI_module.cc new file mode 100644 index 000000000..d8f5062ac --- /dev/null +++ b/sbncode/HitFinder/WireToChannelROI_module.cc @@ -0,0 +1,184 @@ +//////////////////////////////////////////////////////////////////////// +// +// WireToChannelROI class - Convert recob::Wire to recob::ChannelROI +// +// usher@slac.stanford.edu +// +//////////////////////////////////////////////////////////////////////// + +// C/C++ standard libraries +#include +#include +#include // std::pair<> +#include // std::unique_ptr<> + +// framework libraries +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "canvas/Utilities/Exception.h" +#include "canvas/Utilities/InputTag.h" + +// LArSoft libraries +#include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t +#include "larcore/Geometry/WireReadout.h" +#include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom() +#include "larcorealg/CoreUtils/zip.h" +#include "lardataobj/RecoBase/Wire.h" + +#include "sbnobj/ICARUS/TPC/ChannelROI.h" +#include "sbnobj/Common/Utilities/ChannelROICreator.h" + +///creation of calibrated signals on wires +namespace caldata { + +class WireToChannelROI : public art::EDProducer +{ +public: +// create calibrated signals on wires. this class runs +// an fft to remove the electronics shaping. + explicit WireToChannelROI(fhicl::ParameterSet const& pset); + void produce(art::Event& evt); + void beginJob(); + void endJob(); + void reconfigure(fhicl::ParameterSet const& p); +private: + + std::vector fWireModuleLabelVec; ///< vector of modules that made digits + std::vector fOutInstanceLabelVec; ///< The output instance labels to apply + short int fMaxADCScaleFactor; ///< Maximum scale factor to apply + bool fDiagnosticOutput; ///< secret diagnostics flag + size_t fEventCount; ///< count of event processed + + const geo::WireReadoutGeom* fChannelMapAlg = &art::ServiceHandle()->Get(); + +}; // class WireToChannelROI + +DEFINE_ART_MODULE(WireToChannelROI) + +//------------------------------------------------- +WireToChannelROI::WireToChannelROI(fhicl::ParameterSet const& pset) : EDProducer{pset} +{ + this->reconfigure(pset); + + for(const auto& wireLabel : fOutInstanceLabelVec) + { + produces>(wireLabel); + } +} + +////////////////////////////////////////////////////// +void WireToChannelROI::reconfigure(fhicl::ParameterSet const& pset) +{ + // Recover the parameters + fWireModuleLabelVec = pset.get>("WireModuleLabelVec", std::vector()={"decon1droi"}); + fOutInstanceLabelVec = pset.get> ("OutInstanceLabelVec", {"PHYSCRATEDATA"}); + fMaxADCScaleFactor = pset.get("MaxADCScaleFactor", 512); + fDiagnosticOutput = pset.get< bool >("DiagnosticOutput", false); + + if (fWireModuleLabelVec.size() != fOutInstanceLabelVec.size()) + { + throw art::Exception(art::errors::Configuration) << " Configured " << fOutInstanceLabelVec.size() + << " instance names (`OutInstanceLabelVec`) for " << fWireModuleLabelVec.size() + << " input products (`WireModuleLabelVec`)\n"; + } + + return; +} + +//------------------------------------------------- +void WireToChannelROI::beginJob() +{ + fEventCount = 0; +} // beginJob + +////////////////////////////////////////////////////// +void WireToChannelROI::endJob() +{ +} + +////////////////////////////////////////////////////// +void WireToChannelROI::produce(art::Event& evt) +{ + // We need to loop through the list of Wire data we have been given + // This construct from Gianluca Petrillo who invented it and should be given all credit for it! + for(auto const& [wireLabel, instanceName] : util::zip(fWireModuleLabelVec, fOutInstanceLabelVec)) + { + // make a collection of Wires + std::unique_ptr> channelROICol = std::make_unique>(); + + mf::LogInfo("WireToChannelROI") << "WireToChannelROI, looking for Wire data at " << wireLabel.encode(); + + // Read in the collection of full length deconvolved waveforms + const std::vector& wireVec = evt.getProduct>(wireLabel); + + mf::LogInfo("WireToChannelROI") << "--> Recovered Wire data, size: " << wireVec.size(); + + if (!wireVec.empty()) + { + // Reserve the room for the output + channelROICol->reserve(wireVec.size()); + + // Loop through the input ChannelROI collection + for(const auto& wire : wireVec) + { + short int ADCScaleFactor(recob::ChannelROI::defADCScaleFactor); + + // Recover the channel and the view + raw::ChannelID_t channel = wire.Channel(); + + // Create an ROI vector for output + recob::ChannelROI::RegionsOfInterest_t ROIVec; + + // Return the entire waveform (zero padded) so we can find maximum range + std::vector fullWaveform = wire.Signal(); + + std::pair::const_iterator,std::vector::const_iterator> minMaxPair = std::minmax_element(fullWaveform.begin(),fullWaveform.end()); + + float maxValue = std::max(-*minMaxPair.first,*minMaxPair.second); + + ADCScaleFactor = std::max(short(std::round(std::min(float(fMaxADCScaleFactor),static_cast(std::numeric_limits::max())/maxValue))),ADCScaleFactor); + + // Loop through the ROIs for this channel + const recob::Wire::RegionsOfInterest_t& wireROIs = wire.SignalROI(); + + for(const auto& range : wireROIs.get_ranges()) + { + size_t startTick = range.begin_index(); + + std::vector dataVec(range.data().size()); + + for(size_t binIdx = 0; binIdx < range.data().size(); binIdx++) + { + short int scaledIntADCVal = + std::min(std::max(std::round(range.data()[binIdx] * ADCScaleFactor),static_cast(std::numeric_limits::min())), + static_cast(std::numeric_limits::max())); + + dataVec[binIdx] = scaledIntADCVal; + } + + ROIVec.add_range(startTick, std::move(dataVec)); + } + + channelROICol->push_back(recob::ChannelROICreator(std::move(ROIVec),channel,ADCScaleFactor).move()); + } + + mf::LogInfo("WireToChannelROI") << "--> Outputting ChannelROI data, size: " << channelROICol->size() << " with instance name: " << instanceName; + + // Time to stroe everything + if(channelROICol->empty()) mf::LogWarning("WireToChannelROI") << "No wires made for this event."; + } + + evt.put(std::move(channelROICol), instanceName); + } + + fEventCount++; + + return; +} // produce + +} // end namespace caldata diff --git a/sbncode/HitFinder/hitfindermodules_sbn.fcl b/sbncode/HitFinder/hitfindermodules_sbn.fcl new file mode 100644 index 000000000..868e5a2b8 --- /dev/null +++ b/sbncode/HitFinder/hitfindermodules_sbn.fcl @@ -0,0 +1,54 @@ +// from larreco/RecoAlg: +#include "hitalgorithms.fcl" +#include "hitfindermodules.fcl" + +BEGIN_PROLOG + +# This defines a version of the Gauss Hit Finder that uses as input ChannelROIs instead of Wire +# In principle all default parameters remain the same. + +gauss_hitfinder: +{ + module_type: "GaussHitFinderSBN" + CalDataModuleLabel: "caldata" + MaxMultiHit: 5 # maximum hits for multi gaussia fit attempt + AreaMethod: 0 # 0 = area by integral, 1 = area by gaussian area formula + AreaNorms: [ 1.0, 1.0, 1.0 ] # normalizations that put signal area in + # same scale as peak height. + TryNplus1Fits: false # Don't try to refit with extra peak if bad chisq + LongMaxHits: [ 25, 25, 25] # max number hits in long pulse trains + LongPulseWidth: [ 10, 10, 10] # max widths for hits in long pulse trains + Chi2NDF: 500. # maximum Chisquared / NDF allowed to store fit, if fail + # will use "long" pulse method to return hit + AllHitsInstanceName: "" # If non-null then this will be the instance name of all hits output to event + # in this case there will be two hit collections, one filtered and one containing all hits + + # Candididate peak finding done by tool, one tool instantiated per plane (but could be other divisions too) + HitFinderToolVec: + { + CandidateHitsPlane0: @local::candhitfinder_standard # plane 0 + CandidateHitsPlane1: @local::candhitfinder_standard # plane 1 + CandidateHitsPlane2: @local::candhitfinder_standard # plane 2 + } + + # Declare the peak fitting tool + PeakFitter: @local::peakfitter_mrqdt + + # The below are for the hit filtering section of the gaushit finder + FilterHits: false # true = do not keep undesired hits according to settings of HitFilterAlg object + HitFilterAlg: + { + AlgName: "HitFilterAlg" + MinPulseHeight: [5.0, 5.0, 5.0] #minimum hit peak amplitude per plane + MinPulseSigma: [1.0, 1.0, 1.0] #minimum hit rms per plane + } + # In addition to the filter alg we can also filter hits on the same pulse train + PulseHeightCuts: [3.0, 3.0, 3.0 ] # Minimum pulse height + PulseWidthCuts: [2.0, 1.5, 1.0 ] # Minimum pulse width + PulseRatioCuts: [0.35, 0.40, 0.20] # Ratio of pulse height to width +} + +# Define sbn version of gaushit finder +gausshit_sbn: @local::gauss_hitfinder + +END_PROLOG diff --git a/sbncode/HitFinder/wirechannelroiconverters_sbn.fcl b/sbncode/HitFinder/wirechannelroiconverters_sbn.fcl new file mode 100644 index 000000000..a79623b16 --- /dev/null +++ b/sbncode/HitFinder/wirechannelroiconverters_sbn.fcl @@ -0,0 +1,30 @@ +// from larreco/RecoAlg: +#include "hitalgorithms.fcl" +#include "hitfindermodules.fcl" + +BEGIN_PROLOG + +# This defines a version of the Gauss Hit Finder that uses as input ChannelROIs instead of Wire +# In principle all default parameters remain the same. + +channelroitowire: +{ + module_type: "ChannelROIToWire" + WireModuleLabelVec: [" "] + OutInstanceLabelVec: [" "] + DiagnosticOutput: false +} + +wiretochannelroi: +{ + module_type: "WireToChannelROI" + WireModuleLabelVec: [" "] + OutInstanceLabelVec: [" "] + DiagnosticOutput: false +} + +# Define sbn versions +channelroitowire_sbn: @local::channelroitowire +wiretochannelroi_sbn: @local::wiretochannelroi + +END_PROLOG