diff --git a/icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc b/icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc index 438eeaab3..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) { @@ -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 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/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..14bc2e42d 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" @@ -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 @@ -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()->GetProvider()) { this->configure(pset); @@ -228,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.}}); @@ -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 widVec = fChannelMapAlg->ChannelToWire(fChannelIDVec[idx]); @@ -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) @@ -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(), diff --git a/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl b/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl index 525a98204..7a7e2374d 100644 --- a/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl +++ b/icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl @@ -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] ]