Skip to content

Removing bad channels from use in coherent noise corrections #818

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
leoaliaga merged 7 commits into from
May 13, 2025
Merged

Removing bad channels from use in coherent noise corrections #818

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
2016c47
Updating to mark channels listed as bad in the channel status databas…
SFBayLaser Apr 5, 2025
0c58175
Updates aimed at protecting against bad channels
SFBayLaser Apr 23, 2025
80eb1a0
Merge branch 'develop' into feature/usher_markbadchannelsforcoherentn…
SFBayLaser Apr 23, 2025
c069630
Make the output of "ROIs" optional since standard workflows don't use…
SFBayLaser Apr 24, 2025
3403b49
Adding ability to control through fhicl the removal of bad channels (…
SFBayLaser Apr 24, 2025
5cd78cd
Do not label high rms channels as bad
SFBayLaser May 9, 2025
e6dd534
Fix typo - inserting # to make sure comment is a comment
SFBayLaser May 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 24 additions & 10 deletions icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -160,6 +160,7 @@ class DaqDecoderICARUSTPCwROI : public art::ReplicatedProducer

// Fcl parameters.
std::vector<art::InputTag> 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?
Expand Down Expand Up @@ -238,7 +239,9 @@ DaqDecoderICARUSTPCwROI::DaqDecoderICARUSTPCwROI(fhicl::ParameterSet const & pse
for(const auto& fragmentLabel : fFragmentsLabelVec)
{
produces<std::vector<raw::RawDigit>>(fragmentLabel.instance());
produces<std::vector<recob::ChannelROI>>(fragmentLabel.instance());

if (fOutputChannelROIs)
produces<std::vector<recob::ChannelROI>>(fragmentLabel.instance());

if (fOutputRawWaveform)
produces<std::vector<raw::RawDigit>>(fragmentLabel.instance() + fOutputRawWavePath);
Expand Down Expand Up @@ -309,6 +312,7 @@ DaqDecoderICARUSTPCwROI::~DaqDecoderICARUSTPCwROI()
void DaqDecoderICARUSTPCwROI::configure(fhicl::ParameterSet const & pset)
{
fFragmentsLabelVec = pset.get<std::vector<art::InputTag>>("FragmentsLabelVec", std::vector<art::InputTag>()={"daq:PHYSCRATEDATA"});
fOutputChannelROIs = pset.get<bool >("OutputChannelROIs", true);
fOutputRawWaveform = pset.get<bool >("OutputRawWaveform", false);
fOutputCorrection = pset.get<bool >("OutputCorrection", false);
fOutputMorphed = pset.get<bool >("OutputMorphed", false);
Expand Down Expand Up @@ -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::vector<recob::ChannelROI>>(std::move_iterator(concurrentROIs.begin()),
if (fOutputChannelROIs)
{
// Do the same to output the candidate ROIs
ChannelROICollectionPtr channelROICollection = std::make_unique<std::vector<recob::ChannelROI>>(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)
{
Expand Down Expand Up @@ -652,15 +657,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<float,2> meanPos;
Eigen::Vector<float,2> eigenValues;
Eigen::Matrix<float,2,2> 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
Expand Down
22 changes: 20 additions & 2 deletions icaruscode/Decode/DecoderTools/TPCDecoderFilter1D_tool.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -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<lariov::ChannelStatusService const>()->GetProvider())

{
std::cout << "TPCDecoderFilter1D is calling configure method" << std::endl;
this->configure(pset);
Expand Down Expand Up @@ -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;
Expand Down
37 changes: 32 additions & 5 deletions icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -156,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
Expand Down Expand Up @@ -187,9 +191,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<lariov::ChannelStatusService const>()->GetProvider())
{
this->configure(pset);

Expand Down Expand Up @@ -228,6 +236,8 @@ void TPCNoiseFilter1DMC::configure(fhicl::ParameterSet const &pset)
fUseFFTFilter = pset.get<bool >("UseFFTFilter", true);
fLowFreqCorrection = pset.get<bool >("LowFreqCorrection", true);
fDiagnosticOutput = pset.get<bool >("DiagnosticOutput", false);
fRemoveBadChannels = pset.get<bool >("RemoveBadChannels", true);
fRemoveBadRMS = pset.get<bool >("RemoveBadRMS", true);
fFilterModeVec = pset.get<std::vector<std::string>>("FilterModeVec", std::vector<std::string>()={"e","g","d"});

fFFTSigmaValsVec = pset.get<FloatPairVec >("FFTSigmaVals", FloatPairVec()={{1.5,20.}, {1.5,20.}, {2.0,20.}});
Expand Down Expand Up @@ -314,8 +324,26 @@ 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;

// Keep track of the channel
fChannelIDVec[idx] = channelPlaneVec[idx].first;
fChannelIDVec[idx] = channelID;

// Is this a valid channel and what is its status?
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<geo::WireID> widVec = fChannelMapAlg->ChannelToWire(fChannelIDVec[idx]);
Expand Down Expand Up @@ -367,6 +395,8 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&,
fNumTruncBins[idx],
fRangeBins[idx]);

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);
if (fUseFFTFilter)
Expand Down Expand Up @@ -396,9 +426,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(),
Expand Down
4 changes: 3 additions & 1 deletion icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,10 +89,12 @@ TPCNoiseFilter1DTool: {
NSigmaForTrucation: 3.5
StructuringElement: 16
FilterWindow: 10
Threshold: [1.0,1.0,1.0]
Threshold: [20.0,12.,12.] #--> for PCA: [5.0,3.5,3.5]
FilterModeVec: ["e","g","d"]
LowFreqCorrection: true
UseFFTFilter: false
RemoveBadChannels: 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] ]
Expand Down