From 2016c47527fb2cfc6a4b5445be0e3b2940763c1b Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Sat, 5 Apr 2025 07:51:52 -0700 Subject: [PATCH 1/6] Updating to mark channels listed as bad in the channel status database so they are not used in the computation of the coherent noise factors --- .../DecoderTools/TPCDecoderFilter1D_tool.cc | 22 ++++++++++++++-- .../DecoderTools/TPCNoiseFilter1D_tool.cc | 25 +++++++++++++++++-- 2 files changed, 43 insertions(+), 4 deletions(-) diff --git a/icaruscode/Decode/DecoderTools/TPCDecoderFilter1D_tool.cc b/icaruscode/Decode/DecoderTools/TPCDecoderFilter1D_tool.cc index 9cd94a71c..c4be4f88a 100644 --- a/icaruscode/Decode/DecoderTools/TPCDecoderFilter1D_tool.cc +++ b/icaruscode/Decode/DecoderTools/TPCDecoderFilter1D_tool.cc @@ -19,6 +19,8 @@ // LArSoft includes #include "larcore/Geometry/WireReadout.h" +#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h" +#include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h" #include "sbndaq-artdaq-core/Overlays/ICARUS/PhysCrateFragment.hh" @@ -187,9 +189,13 @@ class TPCDecoderFilter1D : virtual public IDecoderFilter // Keep track of the FFT icarus_signal_processing::FFTFilterFunctionVec fFFTFilterFunctionVec; + // channel status DB + const lariov::ChannelStatusProvider* fChannelStatus; }; -TPCDecoderFilter1D::TPCDecoderFilter1D(fhicl::ParameterSet const &pset) +TPCDecoderFilter1D::TPCDecoderFilter1D(fhicl::ParameterSet const &pset) : + fChannelStatus(&art::ServiceHandle()->GetProvider()) + { std::cout << "TPCDecoderFilter1D is calling configure method" << std::endl; this->configure(pset); @@ -408,8 +414,20 @@ void TPCDecoderFilter1D::process_fragment(detinfo::DetectorClocksData const&, icarus_signal_processing::VectorFloat& pedCorDataVec = fPedCorWaveforms[channelOnBoard]; + // Recover the channel ID + int channelID = channelPlanePairVec[chanIdx].first; + // Keep track of the channel - fChannelIDVec[channelOnBoard] = channelPlanePairVec[chanIdx].first; + fChannelIDVec[channelOnBoard] = channelID; + + // Is this a valid channel and what is its status? + if (fChannelStatus->IsPresent(channelID)) + { + // If the channel is bad then we "protect" the entire channel (it will not be used in noise removal) + // Note that the array has already been cleared before calling this function so no need to set opposite case + if (fChannelStatus->IsBad(channelID)) + std::fill(fSelectVals[channelOnBoard].begin(),fSelectVals[channelOnBoard].end(),true); + } // Handle the filter function to use for this channel unsigned int plane = channelPlanePairVec[chanIdx].second; diff --git a/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc b/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc index 64e57e5fe..3de4ded18 100644 --- a/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc +++ b/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc @@ -19,6 +19,8 @@ // LArSoft includes #include "larcore/Geometry/WireReadout.h" +#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h" +#include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h" #include "sbndaq-artdaq-core/Overlays/ICARUS/PhysCrateFragment.hh" @@ -187,9 +189,13 @@ class TPCNoiseFilter1DMC : virtual public INoiseFilter // Keep track of the FFT icarus_signal_processing::FFTFilterFunctionVec fFFTFilterFunctionVec; + + // channel status DB + const lariov::ChannelStatusProvider* fChannelStatus; }; -TPCNoiseFilter1DMC::TPCNoiseFilter1DMC(fhicl::ParameterSet const &pset) +TPCNoiseFilter1DMC::TPCNoiseFilter1DMC(fhicl::ParameterSet const &pset) : + fChannelStatus(&art::ServiceHandle()->GetProvider()) { this->configure(pset); @@ -314,8 +320,23 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&, icarus_signal_processing::VectorFloat& rawDataVec = fRawWaveforms[idx]; icarus_signal_processing::VectorFloat& pedCorDataVec = fPedCorWaveforms[idx]; + // Recover the channel ID + int channelID = channelPlaneVec[idx].first; + // Keep track of the channel - fChannelIDVec[idx] = channelPlaneVec[idx].first; + fChannelIDVec[idx] = channelID; + + // Is this a valid channel and what is its status? + if (fChannelStatus->IsPresent(channelID)) + { + // If the channel is bad then we "protect" the entire channel (it will not be used in noise removal) + // Note that the array has already been cleared before calling this function so no need to set opposite case + if (fChannelStatus->IsBad(channelID)) + { +// std::cout << "--> Channel:" << channelID << " is marked as bad by the channel status service" << std::endl; + std::fill(fSelectVals[idx].begin(),fSelectVals[idx].end(),true); + } + } // We need to recover info on which plane we have std::vector widVec = fChannelMapAlg->ChannelToWire(fChannelIDVec[idx]); From 0c58175910d0f31fb0abaa3c53477edf2837472e Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Wed, 23 Apr 2025 11:11:36 -0700 Subject: [PATCH 2/6] Updates aimed at protecting against bad channels --- .../Decode/DaqDecoderICARUSTPCwROI_module.cc | 15 ++++++++-- .../DecoderTools/TPCNoiseFilter1D_tool.cc | 28 ++++++++++--------- .../DecoderTools/decoderTools_icarus.fcl | 2 +- 3 files changed, 28 insertions(+), 17 deletions(-) diff --git a/icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc b/icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc index 438eeaab3..c8994aae5 100644 --- a/icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc +++ b/icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc @@ -652,15 +652,24 @@ void DaqDecoderICARUSTPCwROI::processSingleFragment(size_t if (fOutputMorphed) { - const icarus_signal_processing::VectorFloat& corrections = decoderTool->getMorphedWaveforms()[chanIdx]; + const icarus_signal_processing::VectorFloat& morphWaveform = decoderTool->getMorphedWaveforms()[chanIdx]; // Need to convert from float to short int - std::transform(corrections.begin(),corrections.end(),wvfm.begin(),[](const auto& val){return short(std::round(val));}); + std::transform(morphWaveform.begin(),morphWaveform.end(),wvfm.begin(),[](const auto& val){return short(std::round(val));}); + + // Get the morphological waveform mean position and spread from PCA + Eigen::Vector meanPos; + Eigen::Vector eigenValues; + Eigen::Matrix eigenVectors; + + icarus_signal_processing::VectorFloat localCopy = morphWaveform; + + waveformTools.principalComponents(localCopy, meanPos, eigenVectors, eigenValues, 6., 3., false); //ConcurrentRawDigitCol::iterator newRawObjItr = coherentRawDigitCol.emplace_back(channel,wvfm.size(),wvfm); ConcurrentRawDigitCol::iterator newRawObjItr = morphedRawDigitCol.push_back(raw::RawDigit(channel,wvfm.size(),wvfm)); - newRawObjItr->SetPedestal(0.,0.); + newRawObjItr->SetPedestal(meanPos[1],std::sqrt(eigenValues[0])); } // Now determine the pedestal and correct for it diff --git a/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc b/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc index 3de4ded18..07402a26b 100644 --- a/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc +++ b/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc @@ -320,6 +320,9 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&, icarus_signal_processing::VectorFloat& rawDataVec = fRawWaveforms[idx]; icarus_signal_processing::VectorFloat& pedCorDataVec = fPedCorWaveforms[idx]; + // Make sure our selection and ROI arrays are initialized + std::fill(fSelectVals[idx].begin(),fSelectVals[idx].end(),false); + // Recover the channel ID int channelID = channelPlaneVec[idx].first; @@ -327,16 +330,16 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&, fChannelIDVec[idx] = channelID; // Is this a valid channel and what is its status? - if (fChannelStatus->IsPresent(channelID)) - { - // If the channel is bad then we "protect" the entire channel (it will not be used in noise removal) - // Note that the array has already been cleared before calling this function so no need to set opposite case - if (fChannelStatus->IsBad(channelID)) - { -// std::cout << "--> Channel:" << channelID << " is marked as bad by the channel status service" << std::endl; - std::fill(fSelectVals[idx].begin(),fSelectVals[idx].end(),true); - } - } + //if (fChannelStatus->IsPresent(channelID)) + //{ + // // If the channel is bad then we "protect" the entire channel (it will not be used in noise removal) + // // Note that the array has already been cleared before calling this function so no need to set opposite case + // if (fChannelStatus->IsBad(channelID)) + // { +// // std::cout << "--> Channel:" << channelID << " is marked as bad by the channel status service" << std::endl; + // std::fill(fSelectVals[idx].begin(),fSelectVals[idx].end(),true); + // } + //} // We need to recover info on which plane we have std::vector widVec = fChannelMapAlg->ChannelToWire(fChannelIDVec[idx]); @@ -388,6 +391,8 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&, fNumTruncBins[idx], fRangeBins[idx]); + if (fFullRMSVals[idx] < 2.5 || fFullRMSVals[idx] > 40.) std::fill(fSelectVals[idx].begin(),fSelectVals[idx].end(),true); + // Convolve with a filter function //if (fUseFFTFilter) (*fFFTFilterFunctionVec[plane])(pedCorDataVec); if (fUseFFTFilter) @@ -417,9 +422,6 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&, // rawDataVec[ rawDataVec.size()/2] = 16. * std::sqrt(eigenValues(0)); // rawDataVec[1+rawDataVec.size()/2] = -16. * std::sqrt(eigenValues(0)); // } - - // Make sure our selection and ROI arrays are initialized - std::fill(fSelectVals[idx].begin(),fSelectVals[idx].end(),false); } (*denoiser)(fWaveLessCoherent.begin(), diff --git a/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl b/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl index 525a98204..e6584d674 100644 --- a/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl +++ b/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl @@ -89,7 +89,7 @@ TPCNoiseFilter1DTool: { NSigmaForTrucation: 3.5 StructuringElement: 16 FilterWindow: 10 - Threshold: [1.0,1.0,1.0] + Threshold: [5.0,3.5,3.5] #[20.0,12.,12.] FilterModeVec: ["e","g","d"] LowFreqCorrection: true UseFFTFilter: false From c069630649dee82a1845b6975369bdc119a1815d Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Thu, 24 Apr 2025 11:34:15 -0700 Subject: [PATCH 3/6] Make the output of "ROIs" optional since standard workflows don't use them at this time so let's save output space --- .../Decode/DaqDecoderICARUSTPCwROI_module.cc | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc b/icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc index c8994aae5..ffd08ab21 100644 --- a/icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc +++ b/icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc @@ -160,6 +160,7 @@ class DaqDecoderICARUSTPCwROI : public art::ReplicatedProducer // Fcl parameters. std::vector fFragmentsLabelVec; ///< The input artdaq fragment label vector (for more than one) + bool fOutputChannelROIs; ///< Should we output the ROIs we found? bool fOutputRawWaveform; ///< Should we output pedestal corrected (not noise filtered)? bool fOutputCorrection; ///< Should we output the coherent noise correction vectors? bool fOutputMorphed; ///< Should we output the morphological filter vectors? @@ -238,7 +239,9 @@ DaqDecoderICARUSTPCwROI::DaqDecoderICARUSTPCwROI(fhicl::ParameterSet const & pse for(const auto& fragmentLabel : fFragmentsLabelVec) { produces>(fragmentLabel.instance()); - produces>(fragmentLabel.instance()); + + if (fOutputChannelROIs) + produces>(fragmentLabel.instance()); if (fOutputRawWaveform) produces>(fragmentLabel.instance() + fOutputRawWavePath); @@ -309,6 +312,7 @@ DaqDecoderICARUSTPCwROI::~DaqDecoderICARUSTPCwROI() void DaqDecoderICARUSTPCwROI::configure(fhicl::ParameterSet const & pset) { fFragmentsLabelVec = pset.get>("FragmentsLabelVec", std::vector()={"daq:PHYSCRATEDATA"}); + fOutputChannelROIs = pset.get("OutputChannelROIs", true); fOutputRawWaveform = pset.get("OutputRawWaveform", false); fOutputCorrection = pset.get("OutputCorrection", false); fOutputMorphed = pset.get("OutputMorphed", false); @@ -416,14 +420,15 @@ void DaqDecoderICARUSTPCwROI::produce(art::Event & event, art::ProcessingFrame c // Now transfer ownership to the event store event.put(std::move(rawDigitCollection), fragmentLabel.instance()); - // Do the same to output the candidate ROIs - ChannelROICollectionPtr channelROICollection = std::make_unique>(std::move_iterator(concurrentROIs.begin()), + if (fOutputChannelROIs) + { + // Do the same to output the candidate ROIs + ChannelROICollectionPtr channelROICollection = std::make_unique>(std::move_iterator(concurrentROIs.begin()), std::move_iterator(concurrentROIs.end())); + std::sort(channelROICollection->begin(),channelROICollection->end(),[](const auto& left, const auto& right){return left.Channel() < right.Channel();}); - std::sort(channelROICollection->begin(),channelROICollection->end(),[](const auto& left, const auto& right){return left.Channel() < right.Channel();}); - - event.put(std::move(channelROICollection), fragmentLabel.instance()); - + event.put(std::move(channelROICollection), fragmentLabel.instance()); + } if (fOutputRawWaveform) { From 3403b49dac7390d7c3d71c9c8989d9e0a9ec4830 Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Thu, 24 Apr 2025 11:35:37 -0700 Subject: [PATCH 4/6] Adding ability to control through fhicl the removal of bad channels (from the bad channel database) and the removal of channels with "bad" truncated RMS values (e.g. too low to be connected and/or too high to be usable). --- .../DecoderTools/TPCNoiseFilter1D_tool.cc | 26 +++++++++++-------- .../DecoderTools/decoderTools_icarus.fcl | 2 ++ 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc b/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc index 07402a26b..14bc2e42d 100644 --- a/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc +++ b/icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc @@ -158,6 +158,8 @@ class TPCNoiseFilter1DMC : virtual public INoiseFilter bool fUseFFTFilter; //< Turn on/off the use of the FFT filter bool fLowFreqCorrection; //< Apply low frequency noise correction bool fDiagnosticOutput; //< If true will spew endless messages to output + bool fRemoveBadChannels; //< Will remove bad channels from calculation + bool fRemoveBadRMS; //< Will ignore channels with too low or too high rms FloatPairVec fFFTSigmaValsVec; //< Gives the sigmas for the filter function FloatPairVec fFFTCutoffValsVec; //< Gives the cutoffs for the filter function std::string fDenoiserType; //< Describes the specific denoiser to use @@ -234,6 +236,8 @@ void TPCNoiseFilter1DMC::configure(fhicl::ParameterSet const &pset) fUseFFTFilter = pset.get("UseFFTFilter", true); fLowFreqCorrection = pset.get("LowFreqCorrection", true); fDiagnosticOutput = pset.get("DiagnosticOutput", false); + fRemoveBadChannels = pset.get("RemoveBadChannels", true); + fRemoveBadRMS = pset.get("RemoveBadRMS", true); fFilterModeVec = pset.get>("FilterModeVec", std::vector()={"e","g","d"}); fFFTSigmaValsVec = pset.get("FFTSigmaVals", FloatPairVec()={{1.5,20.}, {1.5,20.}, {2.0,20.}}); @@ -330,16 +334,16 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&, fChannelIDVec[idx] = channelID; // Is this a valid channel and what is its status? - //if (fChannelStatus->IsPresent(channelID)) - //{ - // // If the channel is bad then we "protect" the entire channel (it will not be used in noise removal) - // // Note that the array has already been cleared before calling this function so no need to set opposite case - // if (fChannelStatus->IsBad(channelID)) - // { -// // std::cout << "--> Channel:" << channelID << " is marked as bad by the channel status service" << std::endl; - // std::fill(fSelectVals[idx].begin(),fSelectVals[idx].end(),true); - // } - //} + if (fRemoveBadChannels && fChannelStatus->IsPresent(channelID)) + { + // If the channel is bad then we "protect" the entire channel (it will not be used in noise removal) + // Note that the array has already been cleared before calling this function so no need to set opposite case + if (fChannelStatus->IsBad(channelID)) + { +// std::cout << "--> Channel:" << channelID << " is marked as bad by the channel status service" << std::endl; + std::fill(fSelectVals[idx].begin(),fSelectVals[idx].end(),true); + } + } // We need to recover info on which plane we have std::vector widVec = fChannelMapAlg->ChannelToWire(fChannelIDVec[idx]); @@ -391,7 +395,7 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&, fNumTruncBins[idx], fRangeBins[idx]); - if (fFullRMSVals[idx] < 2.5 || fFullRMSVals[idx] > 40.) std::fill(fSelectVals[idx].begin(),fSelectVals[idx].end(),true); + if (fRemoveBadRMS && (fFullRMSVals[idx] < 2.5 || fFullRMSVals[idx] > 40.)) std::fill(fSelectVals[idx].begin(),fSelectVals[idx].end(),true); // Convolve with a filter function //if (fUseFFTFilter) (*fFFTFilterFunctionVec[plane])(pedCorDataVec); diff --git a/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl b/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl index e6584d674..a79a1a974 100644 --- a/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl +++ b/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl @@ -93,6 +93,8 @@ TPCNoiseFilter1DTool: { FilterModeVec: ["e","g","d"] LowFreqCorrection: true UseFFTFilter: false + RemoveBadChannels: true + RemoveBadRMS: true FFTSigmaVals: [ [1.5,20.], [1.5,20.], [1.5,20.] ] FFTCutoffVals: [ [8.,800.], [8.,800.], [4.0,800.] ] FragmentIDMap: [ [0,0x140C], [1,0x140E], [2,0x1410], [6,0x1414], [8,0x150E], [9,0x1510] ] From 5cd78cd6e271b50c2f395f44b61c170712e2655f Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Fri, 9 May 2025 10:44:46 -0700 Subject: [PATCH 5/6] Do not label high rms channels as bad --- icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl b/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl index a79a1a974..5a25b690f 100644 --- a/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl +++ b/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl @@ -89,12 +89,12 @@ TPCNoiseFilter1DTool: { NSigmaForTrucation: 3.5 StructuringElement: 16 FilterWindow: 10 - Threshold: [5.0,3.5,3.5] #[20.0,12.,12.] + Threshold: [20.0,12.,12.] --> for PCA: [5.0,3.5,3.5] FilterModeVec: ["e","g","d"] LowFreqCorrection: true UseFFTFilter: false RemoveBadChannels: true - RemoveBadRMS: true + RemoveBadRMS: false FFTSigmaVals: [ [1.5,20.], [1.5,20.], [1.5,20.] ] FFTCutoffVals: [ [8.,800.], [8.,800.], [4.0,800.] ] FragmentIDMap: [ [0,0x140C], [1,0x140E], [2,0x1410], [6,0x1414], [8,0x150E], [9,0x1510] ] From e6dd5348faecb32856ffbed0c3964112e4649dc4 Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Fri, 9 May 2025 11:32:54 -0700 Subject: [PATCH 6/6] Fix typo - inserting # to make sure comment is a comment --- icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl b/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl index 5a25b690f..7a7e2374d 100644 --- a/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl +++ b/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl @@ -89,7 +89,7 @@ TPCNoiseFilter1DTool: { NSigmaForTrucation: 3.5 StructuringElement: 16 FilterWindow: 10 - Threshold: [20.0,12.,12.] --> for PCA: [5.0,3.5,3.5] + Threshold: [20.0,12.,12.] #--> for PCA: [5.0,3.5,3.5] FilterModeVec: ["e","g","d"] LowFreqCorrection: true UseFFTFilter: false