From ca7e8344fbd193e7ecdc0afd241430b57c2c42b8 Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Tue, 22 Apr 2025 13:31:51 -0700 Subject: [PATCH] The gauss hit finder operating on ChannelROIs has moved to sbncode, as well the modules converting from Wire to/from ChannelROI. Removing from icaruscode --- fcl/reco/Definitions/stage0_icarus_defs.fcl | 9 +- .../SignalProcessing/HitFinder/CMakeLists.txt | 6 - .../HitFinder/GaussHitFinderICARUS_module.cc | 628 ------------------ .../HitFinder/hitfindermodules_icarus.fcl | 39 +- .../SignalProcessing/RecoWire/CMakeLists.txt | 1 - .../RecoWire/ROIConverter_module.cc | 165 ----- 6 files changed, 19 insertions(+), 829 deletions(-) delete mode 100644 icaruscode/TPC/SignalProcessing/HitFinder/GaussHitFinderICARUS_module.cc delete mode 100644 icaruscode/TPC/SignalProcessing/RecoWire/ROIConverter_module.cc diff --git a/fcl/reco/Definitions/stage0_icarus_defs.fcl b/fcl/reco/Definitions/stage0_icarus_defs.fcl index c22d2f9a9..e832aa6ba 100644 --- a/fcl/reco/Definitions/stage0_icarus_defs.fcl +++ b/fcl/reco/Definitions/stage0_icarus_defs.fcl @@ -7,6 +7,7 @@ #include "decoderdefs_icarus.fcl" #include "recowire_icarus.fcl" +#include "wirechannelroiconverters_sbn.fcl" #include "hitfindermodules_icarus.fcl" #include "timing_icarus.fcl" #include "timing_beam.fcl" @@ -90,10 +91,10 @@ icarus_stage0_producers: gaushit2dTPCEW: @local::gaus_hitfinder_icarus gaushit2dTPCEE: @local::gaus_hitfinder_icarus - gausshitTPCWW: @local::gauss_hitfinder_icarus - gausshitTPCWE: @local::gauss_hitfinder_icarus - gausshitTPCEW: @local::gauss_hitfinder_icarus - gausshitTPCEE: @local::gauss_hitfinder_icarus + gausshitTPCWW: @local::gausshit_sbn + gausshitTPCWE: @local::gausshit_sbn + gausshitTPCEW: @local::gausshit_sbn + gausshitTPCEE: @local::gausshit_sbn ### CRT hit finder producer crthit: @local::standard_crthitproducer diff --git a/icaruscode/TPC/SignalProcessing/HitFinder/CMakeLists.txt b/icaruscode/TPC/SignalProcessing/HitFinder/CMakeLists.txt index 17d747916..06e6714aa 100644 --- a/icaruscode/TPC/SignalProcessing/HitFinder/CMakeLists.txt +++ b/icaruscode/TPC/SignalProcessing/HitFinder/CMakeLists.txt @@ -46,12 +46,6 @@ cet_build_plugin(HitMerger art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(HitSelector art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(ICARUSHitFinder art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(HitConverter art::module LIBRARIES ${MODULE_LIBRARIES}) -cet_build_plugin(GaussHitFinderICARUS art::module - LIBRARIES ${MODULE_LIBRARIES} - larreco::HitFinder - larreco::CandidateHitFinderTool - larreco::PeakFitterTool - ) install_headers() diff --git a/icaruscode/TPC/SignalProcessing/HitFinder/GaussHitFinderICARUS_module.cc b/icaruscode/TPC/SignalProcessing/HitFinder/GaussHitFinderICARUS_module.cc deleted file mode 100644 index 240041172..000000000 --- a/icaruscode/TPC/SignalProcessing/HitFinder/GaussHitFinderICARUS_module.cc +++ /dev/null @@ -1,628 +0,0 @@ -//////////////////////////////////////////////////////////////////////// -// -// GaussHitFinderICARUS 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_GaussHitFinder -// gaushit: @local::argoneut_GaussHitFinder -//////////////////////////////////////////////////////////////////////// - -// 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 "lardataobj/RecoBase/Hit.h" -#include "larreco/HitFinder/HitFilterAlg.h" -#include "icaruscode/TPC/Utilities/HitCreator.h" - -#include "sbnobj/ICARUS/TPC/ChannelROI.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 GaussHitFinderICARUS : public art::SharedProducer { - public: - explicit GaussHitFinderICARUS(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 GaussHitFinderICARUS - - //------------------------------------------------- - //------------------------------------------------- - GaussHitFinderICARUS::GaussHitFinderICARUS(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) - icarus::HitCollectionCreator::declare_products( - producesCollector(), fAllHitsInstanceName, true, false); //fMakeRawDigitAssns); - - // and now the filtered hits... - if (fAllHitsInstanceName != "") - icarus::HitCollectionCreator::declare_products( - producesCollector(), "", true, false); //fMakeRawDigitAssns); - - return; - } // GaussHitFinderICARUS::GaussHitFinderICARUS() - - //------------------------------------------------- - //------------------------------------------------- - std::vector GaussHitFinderICARUS::FillOutHitParameterVector(const std::vector& input) - { - if (input.size() == 0) - throw std::runtime_error( - "GaussHitFinderICARUS::FillOutHitParameterVector ERROR! Input config vector has zero size."); - - std::vector output; - const unsigned int N_PLANES = art::ServiceHandle()->Get().Nplanes(); - - if (input.size() == 1) - output.resize(N_PLANES, input[0]); - else if (input.size() == N_PLANES) - output = input; - else - throw std::runtime_error("GaussHitFinderICARUS::FillOutHitParameterVector ERROR! Input config " - "vector size !=1 and !=N_PLANES."); - return output; - } - - //------------------------------------------------- - //------------------------------------------------- - void GaussHitFinderICARUS::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 GaussHitFinderICARUS::produce(art::Event& evt, art::ProcessingFrame const&) - { - unsigned int count = fEventCount.fetch_add(1); - //================================================================================================== - - TH1::AddDirectory(kFALSE); - - // Instantiate and Reset a stop watch - //TStopwatch StopWatch; - //StopWatch.Reset(); - - // ################################ - // ### Calling Geometry service ### - // ################################ - auto const& wireReadout = 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 - icarus::HitCollectionCreator allHitCol(evt, fAllHitsInstanceName, true, false); - - // Handle the filtered hits collection... - icarus::HitCollectionCreator hcol(evt, "", true, false); - icarus::HitCollectionCreator* filteredHitCol = 0; - - if (fFilterHits) filteredHitCol = &hcol; - - //store in a thread safe way - struct hitstruct { - icarus::Hit hit_tbb; - art::Ptr wire_tbb; - }; - - tbb::concurrent_vector hitstruct_vec; - tbb::concurrent_vector filthitstruct_vec; - - // if (fAllHitsInstanceName != "") filteredHitCol = &hcol; - - // ########################################## - // ### Reading in the Wire List object(s) ### - // ########################################## - art::Handle> wireVecHandle; - evt.getByLabel(fCalDataModuleLabel, wireVecHandle); - - //################################################# - //### 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; - }; - - - std::cout << " " << std::endl; - std::cout << "************************** In gauss hit finder ***************************" << std::endl; - - //############################## - //### Looping over the wires ### - //############################## - //for(size_t wireIter = 0; wireIter < wireVecHandle->size(); wireIter++) - //{ - tbb::parallel_for( - static_cast(0), - wireVecHandle->size(), - [&](size_t& wireIter) { - // #################################### - // ### Getting this particular wire ### - // #################################### - art::Ptr wire(wireVecHandle, wireIter); - - // --- Setting Channel Number and Signal type --- - - raw::ChannelID_t channel = wire->Channel(); - - // get the WireID for this hit - std::vector wids = wireReadout.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_t& signalROI = wire->SignalROI(); - - // This will be a handy definition - using SparseVectorFloat_t = lar::sparse_vector; - - // for (const auto& range : signalROI.get_ranges()) { - tbb::parallel_for( - static_cast(0), - signalROI.n_ranges(), - [&](size_t& rangeIter) { - const auto& rangeShort = signalROI.range(rangeIter); - // ROI start time - raw::TDCtick_t roiFirstBinTick = rangeShort.begin_index(); - - // For testing we are just going to copy to a float sparse vector and go from there - std::vector floatADCvec(rangeShort.data().size()); - - std::transform(rangeShort.data().begin(),rangeShort.data().end(),floatADCvec.begin(),[](short adc){return static_cast(adc);}); - - SparseVectorFloat_t tempSparseVector; - - tempSparseVector.add_range(rangeShort.begin_index(),floatADCvec); - - const auto& range = tempSparseVector.get_ranges().front(); - - // ########################################################### - // ### 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 sumADC = - std::accumulate(range.begin() + firstTick, range.begin() + lastTick, 0.); - double peakSigma = (lastTick - firstTick) / 3.; // Set the width... - double peakAmp = 0.3989 * sumADC / 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; - - for (const auto& peakParams : peakParamsVec) { - // Extract values for this hit - float peakAmp = peakParams.peakAmplitude; - float peakMean = peakParams.peakCenter; - float peakWidth = peakParams.peakSigma; - - // Place one bit of protection here - if (std::isnan(peakAmp)) { - std::cout << "**** hit peak amplitude is a nan! Channel: " << channel - << ", start tick: " << startT << std::endl; - continue; - } - - // Extract errors - float peakAmpErr = peakParams.peakAmplitudeError; - float peakMeanErr = peakParams.peakCenterError; - float peakWidthErr = peakParams.peakSigmaError; - - // ### Charge ### - float charge = - chargeFunc(peakMean, peakAmp, peakWidth, fAreaNormsVec[plane], startT, endT); - ; - float chargeErr = - std::sqrt(TMath::Pi()) * (peakAmpErr * peakWidthErr + peakWidthErr * peakAmpErr); - - // ### limits for getting sums - std::vector::const_iterator sumStartItr = range.begin() + startT; - std::vector::const_iterator sumEndItr = range.begin() + endT; - - // ### Sum of ADC counts - double sumADC = std::accumulate(sumStartItr, sumEndItr, 0.); - - // ok, now create the hit - icarus::HitCreator hitcreator( - *wire, // wire reference - wid, // wire ID - startT + roiFirstBinTick, // start_tick TODO check - endT + roiFirstBinTick, // end_tick TODO check - peakWidth, // rms - peakMean + roiFirstBinTick, // peak_time - peakMeanErr, // sigma_peak_time - peakAmp, // peak_amplitude - peakAmpErr, // sigma_peak_amplitude - charge, // hit_integral - chargeErr, // hit_sigma_integral - sumADC, // summedADC FIXME - nGausForFit, // multiplicity - numHits, // local_index TODO check that the order is correct - chi2PerNDF, // goodness_of_fit - NDF // dof - ); - - if (filteredHitCol) filteredHitVec.push_back(hitcreator.copy()); - - const icarus::Hit hit(hitcreator.move()); - - // This loop will store ALL hits - hitstruct tmp{std::move(hit), wire}; - 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), wire}; -// 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].wire_tbb); - } - - for (size_t j = 0; j < filthitstruct_vec.size(); j++) { - filteredHitCol->emplace_back(filthitstruct_vec[j].hit_tbb, filthitstruct_vec[j].wire_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); - } - - // Keep track of events processed - //fEventCount++; - - } // End of produce() - - DEFINE_ART_MODULE(GaussHitFinderICARUS) - -} // end of hit namespace diff --git a/icaruscode/TPC/SignalProcessing/HitFinder/hitfindermodules_icarus.fcl b/icaruscode/TPC/SignalProcessing/HitFinder/hitfindermodules_icarus.fcl index 722f6355c..5df56b46b 100644 --- a/icaruscode/TPC/SignalProcessing/HitFinder/hitfindermodules_icarus.fcl +++ b/icaruscode/TPC/SignalProcessing/HitFinder/hitfindermodules_icarus.fcl @@ -2,6 +2,7 @@ #include "hitalgorithms.fcl" #include "HitFinderTools_ICARUS.fcl" #include "hitfindermodules.fcl" +#include "hitfindermodules_sbn.fcl" BEGIN_PROLOG @@ -100,31 +101,19 @@ gaus_hitfinder_icarus.HitFinderToolVec.CandidateHitsPlane2: gaus_hitfinder_icarus.HitFinderToolVec.CandidateHitsPlane2.Plane: 2 gaus_hitfinder_icarus.HitFinderToolVec.CandidateHitsPlane2.RoiThreshold: 9. - -# Define icarus version of test version gausshit finder -gauss_hitfinder_icarus: @local::gaus_hitfinder - -gauss_hitfinder_icarus.module_type: GaussHitFinderICARUS -gauss_hitfinder_icarus.CalDataModuleLabel: "decon1droi" -gauss_hitfinder_icarus.AreaNorms: [ 1.0, 1.0, 1.0 ] -gauss_hitfinder_icarus.MaxMultiHit: 5 -gauss_hitfinder_icarus.TryNplus1Fits: false -gauss_hitfinder_icarus.Chi2NDF: 500. -gauss_hitfinder_icarus.PeakFitter.MinWidth: 1 -gauss_hitfinder_icarus.PeakFitter.FloatBaseline: false -gauss_hitfinder_icarus.PeakFitter.tool_type: "PeakFitterMrqdt" -gauss_hitfinder_icarus.LongMaxHits: [25, 25, 25] -gauss_hitfinder_icarus.LongPulseWidth: [10, 10, 10] -gauss_hitfinder_icarus.HitFinderToolVec.CandidateHitsPlane0: @local::candhitfinder_standard -gauss_hitfinder_icarus.HitFinderToolVec.CandidateHitsPlane0.Plane: 0 -gauss_hitfinder_icarus.gaushit2dTPCWW.HitFinderToolVec.CandidateHitsPlane0.RoiThreshold: 9. -gauss_hitfinder_icarus.HitFinderToolVec.CandidateHitsPlane1: @local::candhitfinder_standard -gauss_hitfinder_icarus.HitFinderToolVec.CandidateHitsPlane1.Plane: 1 -gauss_hitfinder_icarus.gaushit2dTPCWW.HitFinderToolVec.CandidateHitsPlane0.RoiThreshold: 9.5 -gauss_hitfinder_icarus.HitFinderToolVec.CandidateHitsPlane2: @local::candhitfinder_standard -gauss_hitfinder_icarus.HitFinderToolVec.CandidateHitsPlane2.Plane: 2 -gauss_hitfinder_icarus.gaushit2dTPCWW.HitFinderToolVec.CandidateHitsPlane0.RoiThreshold: 9. - +# Now set parameters for the GaussHitFinderSBN +# Note these are "default" for ICARUS and will move to a fcl file in icaruscode +gausshit_sbn.PeakFitter.MinWidth: 1 +gausshit_sbn.PeakFitter.FloatBaseline: false +gausshit_sbn.HitFinderToolVec.CandidateHitsPlane0: @local::candhitfinder_standard +gausshit_sbn.HitFinderToolVec.CandidateHitsPlane0.Plane: 0 +gausshit_sbn.HitFinderToolVec.CandidateHitsPlane0.RoiThreshold: 9. +gausshit_sbn.HitFinderToolVec.CandidateHitsPlane1: @local::candhitfinder_standard +gausshit_sbn.HitFinderToolVec.CandidateHitsPlane1.Plane: 1 +gausshit_sbn.HitFinderToolVec.CandidateHitsPlane1.RoiThreshold: 9.5 +gausshit_sbn.HitFinderToolVec.CandidateHitsPlane2: @local::candhitfinder_standard +gausshit_sbn.HitFinderToolVec.CandidateHitsPlane2.Plane: 2 +gausshit_sbn.HitFinderToolVec.CandidateHitsPlane2.RoiThreshold: 9. icarus_hitconverter: { diff --git a/icaruscode/TPC/SignalProcessing/RecoWire/CMakeLists.txt b/icaruscode/TPC/SignalProcessing/RecoWire/CMakeLists.txt index f97590cf6..3454125fb 100644 --- a/icaruscode/TPC/SignalProcessing/RecoWire/CMakeLists.txt +++ b/icaruscode/TPC/SignalProcessing/RecoWire/CMakeLists.txt @@ -50,7 +50,6 @@ set( MODULE_LIBRARIES FFTW3::FFTW3 ) cet_build_plugin(Decon1DROI art::module LIBRARIES ${MODULE_LIBRARIES}) -cet_build_plugin(ROIConverter art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(ROIFinder art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(SimChannelROI art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(RecoWireICARUSRaw art::module LIBRARIES ${MODULE_LIBRARIES}) diff --git a/icaruscode/TPC/SignalProcessing/RecoWire/ROIConverter_module.cc b/icaruscode/TPC/SignalProcessing/RecoWire/ROIConverter_module.cc deleted file mode 100644 index c53c17e8a..000000000 --- a/icaruscode/TPC/SignalProcessing/RecoWire/ROIConverter_module.cc +++ /dev/null @@ -1,165 +0,0 @@ -//////////////////////////////////////////////////////////////////////// -// -// ROIConvert class - An ROI finding module for complete deconvolved waveforms -// -// 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 ROIConvert : public art::EDProducer -{ -public: -// create calibrated signals on wires. this class runs -// an fft to remove the electronics shaping. - explicit ROIConvert(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 ROIConvert - -DEFINE_ART_MODULE(ROIConvert) - -//------------------------------------------------- -ROIConvert::ROIConvert(fhicl::ParameterSet const& pset) : EDProducer{pset} -{ - this->reconfigure(pset); - - for(const auto& wireLabel : fOutInstanceLabelVec) - { - produces>(wireLabel); - } -} - -////////////////////////////////////////////////////// -void ROIConvert::reconfigure(fhicl::ParameterSet const& pset) -{ - // Recover the parameters - fWireModuleLabelVec = pset.get>("WireModuleLabelVec", std::vector()={"decon1droi"}); - fOutInstanceLabelVec = pset.get> ("OutInstanceLabelVec", {"PHYSCRATEDATA"}); - fDiagnosticOutput = pset.get< bool >("DaignosticOutput", false); - - if (fWireModuleLabelVec.size() != fWireModuleLabelVec.size()) - { - throw art::Exception(art::errors::Configuration) << " Configured " << fOutInstanceLabelVec.size() - << " instance names (`OutInstanceLabelVec`) for " << fWireModuleLabelVec.size() - << " input products (`WireModuleLabelVec`)\n"; - } - - return; -} - -//------------------------------------------------- -void ROIConvert::beginJob() -{ - fEventCount = 0; -} // beginJob - -////////////////////////////////////////////////////// -void ROIConvert::endJob() -{ -} - -////////////////////////////////////////////////////// -void ROIConvert::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("ROIConvert") << "ROIConvert, 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("ROIConvert") << "--> 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); - - // 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] = range.data()[binIdx]; - - ROIVec.add_range(startTick, std::move(dataVec)); - } - - wireCol->push_back(recob::WireCreator(std::move(ROIVec),channel,view).move()); - } - - mf::LogInfo("ROIConvert") << "--> Outputting Wire data, size: " << wireCol->size() << " with instance name: " << instanceName; - - // Time to stroe everything - if(wireCol->empty()) mf::LogWarning("ROIConvert") << "No wires made for this event."; - } - - evt.put(std::move(wireCol), instanceName); - } - - fEventCount++; - - return; -} // produce - -} // end namespace caldata