Skip to content

Commit

Permalink
Merge pull request #332 from SBNSoftware/develop
Browse files Browse the repository at this point in the history 
Weekly Release
  • Loading branch information
@rennney
rennney authored Feb 4, 2022
2 parents 4383378 + fad699c commit 4403b63
Show file tree
Hide file tree
Showing 15 changed files with 858 additions and 46 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@

cmake_minimum_required(VERSION 3.19 FATAL_ERROR)

project(icaruscode VERSION 09.42.01 LANGUAGES CXX)
project(icaruscode VERSION 09.42.02 LANGUAGES CXX)

message(STATUS
"\n-- ============================================================================="
Expand Down
2 changes: 1 addition & 1 deletion icaruscode/Analysis/TPCPurityInfoAna_module.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ class ana::TPCPurityInfoAna : public art::EDAnalyzer {
ana::TPCPurityInfoAna::TPCPurityInfoAna(fhicl::ParameterSet const& p)
: EDAnalyzer{p} // ,
, fPurityInfoLabel(p.get<art::InputTag>("PurityInfoLabel"))
, fPrintInfo(p.get<bool>("PrintInfo",true))
, fPrintInfo(p.get<bool>("PrintInfo",false))
{
consumes< std::vector<anab::TPCPurityInfo> >(fPurityInfoLabel);
}
Expand Down
708 changes: 708 additions & 0 deletions icaruscode/Analysis/TPCPurityMonitor_module.cc
View file

Large diffs are not rendered by default.

17 changes: 17 additions & 0 deletions icaruscode/Analysis/analysis_icarus.fcl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,5 +33,22 @@ icarus_HitEfficiencyAnalysis:
HitEfficiencyHistogramToolList: [ @local::TrackHitEfficiencyAnalysisTool ]
}

icarus_TPCPurityMonitor:
{
module_type: TPCPurityMonitor
# TrackLabel: ["pandoraKalmanTrackGausCryoW","pandoraKalmanTrackGausCryoE"]
TrackLabel: ["pandoraTrackGausCryoW","pandoraTrackGausCryoE"]
SelectedPlane: 2
MinNumHits: 100
MinTickRange: 150.
AssumedELifetime: 600000.
MinRejectFraction: 0.05
MaxRejectFraction: 0.95
OutlierRejectFrac: 0.70
UseHitIntegral: true
WeightByChiSq: false
DiagnosticTuple: true
}

END_PROLOG

33 changes: 30 additions & 3 deletions icaruscode/Decode/DaqDecoderICARUSTPC_module.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@
#include "icaruscode/Decode/DecoderTools/IDecoderFilter.h"

#include "icarus_signal_processing/ICARUSSigProcDefs.h"
#include "icarus_signal_processing/WaveformTools.h"

namespace daq
{
Expand Down Expand Up @@ -126,6 +127,7 @@ class DaqDecoderICARUSTPC : public art::ReplicatedProducer
std::string fOutputRawWavePath; ///< Path to assign to the output if asked for
std::string fOutputCoherentPath; ///< Path to assign to the output if asked for
unsigned int fPlaneToSimulate; ///< Use to get fragment offset
float fSigmaForTruncation; ///< This determines the point at which we truncate bins for the RMS calc

// Statistics.
int fNumEvent; ///< Number of events seen.
Expand Down Expand Up @@ -234,6 +236,7 @@ void DaqDecoderICARUSTPC::configure(fhicl::ParameterSet const & pset)
fOutputRawWavePath = pset.get<std::string >("OutputRawWavePath", "raw");
fOutputCoherentPath = pset.get<std::string >("OutputCoherentPath", "Cor");
fPlaneToSimulate = pset.get<unsigned int >("PlaneToSimulate", 2);
fSigmaForTruncation = pset.get<float >("NSigmaForTrucation", 3.5);
}

//----------------------------------------------------------------------------
Expand Down Expand Up @@ -376,11 +379,35 @@ void DaqDecoderICARUSTPC::processSingleFragment(size_t

double totalTime = theClockProcess.accumulated_real_time();

// We need to recalculate pedestals for the noise corrected waveforms
icarus_signal_processing::WaveformTools<float> waveformTools;

// Save the filtered RawDigitsactive but for corrected raw digits pedestal is zero
const icarus_signal_processing::VectorFloat locPedsVec(decoderTool->getWaveLessCoherent().size(),0.);
const icarus_signal_processing::VectorInt& channelVec = decoderTool->getChannelIDs();
icarus_signal_processing::VectorFloat locPedsVec(decoderTool->getWaveLessCoherent().size(),0.);
icarus_signal_processing::VectorFloat locFullRMSVec(locPedsVec.size(),0.);
icarus_signal_processing::VectorFloat locTruncRMSVec(locPedsVec.size(),0.);
icarus_signal_processing::VectorInt locNumTruncBins(locPedsVec.size(),0);
icarus_signal_processing::VectorInt locRangeBins(locPedsVec.size(),0);

const icarus_signal_processing::VectorInt& channelVec = decoderTool->getChannelIDs();
const icarus_signal_processing::ArrayFloat& corWaveforms = decoderTool->getWaveLessCoherent();

icarus_signal_processing::ArrayFloat pedCorWaveforms(corWaveforms.size(),icarus_signal_processing::VectorFloat(corWaveforms[0].size()));

for(size_t idx = 0; idx < corWaveforms.size(); idx++)
{
// Now determine the pedestal and correct for it
waveformTools.getPedestalCorrectedWaveform(corWaveforms[idx],
pedCorWaveforms[idx],
fSigmaForTruncation,
locPedsVec[idx],
locFullRMSVec[idx],
locTruncRMSVec[idx],
locNumTruncBins[idx],
locRangeBins[idx]);
}

saveRawDigits(decoderTool->getWaveLessCoherent(),locPedsVec,decoderTool->getTruncRMSVals(), channelVec, rawDigitCollection);
saveRawDigits(pedCorWaveforms, locPedsVec, locTruncRMSVec, channelVec, rawDigitCollection);

// Optionally, save the raw RawDigits
if (fOutputRawWaveform)
Expand Down
52 changes: 40 additions & 12 deletions icaruscode/Decode/DaqDecoderICARUSTPCwROI_module.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,8 @@ class DaqDecoderICARUSTPCwROI : public art::ReplicatedProducer
std::string fOutputRawWavePath; ///< Path to assign to the output if asked for
std::string fOutputCoherentPath; ///< Path to assign to the output if asked for
bool fDiagnosticOutput; ///< Set this to get lots of messages
float fSigmaForTruncation; ///< Cut for truncated rms calc
size_t fCoherentNoiseGrouping; ///< Grouping for removing coherent noise

const std::string fLogCategory; ///< Output category when logging messages

Expand Down Expand Up @@ -292,12 +294,14 @@ DaqDecoderICARUSTPCwROI::~DaqDecoderICARUSTPCwROI()
///
void DaqDecoderICARUSTPCwROI::configure(fhicl::ParameterSet const & pset)
{
fFragmentsLabelVec = pset.get<std::vector<art::InputTag>>("FragmentsLabelVec", std::vector<art::InputTag>()={"daq:PHYSCRATEDATA"});
fOutputRawWaveform = pset.get<bool >("OutputRawWaveform", false);
fOutputCorrection = pset.get<bool >("OutputCorrection", false);
fOutputRawWavePath = pset.get<std::string >("OutputRawWavePath", "raw");
fOutputCoherentPath = pset.get<std::string >("OutputCoherentPath", "Cor");
fDiagnosticOutput = pset.get<bool >("DiagnosticOutput", false);
fFragmentsLabelVec = pset.get<std::vector<art::InputTag>>("FragmentsLabelVec", std::vector<art::InputTag>()={"daq:PHYSCRATEDATA"});
fOutputRawWaveform = pset.get<bool >("OutputRawWaveform", false);
fOutputCorrection = pset.get<bool >("OutputCorrection", false);
fOutputRawWavePath = pset.get<std::string >("OutputRawWavePath", "raw");
fOutputCoherentPath = pset.get<std::string >("OutputCoherentPath", "Cor");
fDiagnosticOutput = pset.get<bool >("DiagnosticOutput", false);
fSigmaForTruncation = pset.get<float >("NSigmaForTrucation", 3.5);
fCoherentNoiseGrouping = pset.get<size_t >("CoherentGrouping", 64);
}

//----------------------------------------------------------------------------
Expand Down Expand Up @@ -571,7 +575,24 @@ void DaqDecoderICARUSTPCwROI::processSingleFragment(size_t
}

//process_fragment(event, rawfrag, product_collection, header_collection);
decoderTool->process_fragment(clockData, channelArrayPair.first, channelArrayPair.second);
decoderTool->process_fragment(clockData, channelArrayPair.first, channelArrayPair.second, fCoherentNoiseGrouping);

// We need to recalculate pedestals for the noise corrected waveforms
icarus_signal_processing::WaveformTools<float> waveformTools;

// Local storage for recomputing the the pedestals for the noise corrected data
float localPedestal(0.);
float localFullRMS(0.);
float localTruncRMS(0.);
int localNumTruncBins(0);
int localRangeBins(0);

float sigmaCut(fSigmaForTruncation);

// Recover the denoised waveform
const icarus_signal_processing::ArrayFloat& denoised = decoderTool->getWaveLessCoherent();

icarus_signal_processing::VectorFloat pedCorWaveforms(denoised[0].size());

for(size_t chanIdx = 0; chanIdx < nChannelsPerBoard; chanIdx++)
{
Expand Down Expand Up @@ -604,15 +625,22 @@ void DaqDecoderICARUSTPCwROI::processSingleFragment(size_t
newRawObjItr->SetPedestal(0.,0.);
}

// Recover the denoised waveform
const icarus_signal_processing::VectorFloat& denoised = decoderTool->getWaveLessCoherent()[chanIdx];
// Now determine the pedestal and correct for it
waveformTools.getPedestalCorrectedWaveform(denoised[chanIdx],
pedCorWaveforms,
sigmaCut,
localPedestal,
localFullRMS,
localTruncRMS,
localNumTruncBins,
localRangeBins);

// Need to convert from float to short int
std::transform(denoised.begin(),denoised.end(),wvfm.begin(),[](const auto& val){return short(std::round(val));});
std::transform(denoised[chanIdx].begin(),denoised[chanIdx].end(),wvfm.begin(),[](const auto& val){return short(std::round(val));});

ConcurrentRawDigitCol::iterator newObjItr = concurrentRawDigitCol.emplace_back(channel,wvfm.size(),wvfm);

newObjItr->SetPedestal(0.,decoderTool->getTruncRMSVals()[chanIdx]);
newObjItr->SetPedestal(localPedestal,localFullRMS);

// And, finally, the ROIs
const icarus_signal_processing::VectorBool& chanROIs = decoderTool->getROIVals()[chanIdx];
Expand All @@ -631,7 +659,7 @@ void DaqDecoderICARUSTPCwROI::processSingleFragment(size_t
{
std::vector<float> holder(roiIdx - roiStartIdx);

for(size_t idx = 0; idx < holder.size(); idx++) holder[idx] = denoised[roiStartIdx+idx];
for(size_t idx = 0; idx < holder.size(); idx++) holder[idx] = denoised[chanIdx][roiStartIdx+idx];

ROIVec.add_range(roiStartIdx, std::move(holder));
}
Expand Down
3 changes: 2 additions & 1 deletion icaruscode/Decode/DecoderTools/INoiseFilter.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,8 @@ class INoiseFilter

virtual void process_fragment(detinfo::DetectorClocksData const&,
const daq::INoiseFilter::ChannelPlaneVec&,
const icarus_signal_processing::ArrayFloat&) = 0;
const icarus_signal_processing::ArrayFloat&,
const size_t&) = 0;

/**
* @brief Recover the channels for the processed fragment
Expand Down
17 changes: 9 additions & 8 deletions icaruscode/Decode/DecoderTools/TPCNoiseFilter1D_tool.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,8 @@ class TPCNoiseFilter1DMC : virtual public INoiseFilter
*/
virtual void process_fragment(detinfo::DetectorClocksData const&,
const daq::INoiseFilter::ChannelPlaneVec&,
const icarus_signal_processing::ArrayFloat&) override;
const icarus_signal_processing::ArrayFloat&,
const size_t&) override;

/**
* @brief Recover the channels for the processed fragment
Expand Down Expand Up @@ -141,7 +142,6 @@ class TPCNoiseFilter1DMC : virtual public INoiseFilter
using FloatPairVec = std::vector<std::pair<float,float>>;

float fSigmaForTruncation; //< Selection cut for truncated rms calculation
size_t fCoherentNoiseGrouping; //< # channels in common for coherent noise
size_t fCoherentNoiseOffset; //< Offset for midplane
size_t fStructuringElement; //< Structuring element for morphological filter
size_t fMorphWindow; //< Window for filter
Expand Down Expand Up @@ -212,7 +212,6 @@ TPCNoiseFilter1DMC::~TPCNoiseFilter1DMC()
void TPCNoiseFilter1DMC::configure(fhicl::ParameterSet const &pset)
{
fSigmaForTruncation = pset.get<float >("NSigmaForTrucation", 3.5);
fCoherentNoiseGrouping = pset.get<size_t >("CoherentGrouping", 64);
fCoherentNoiseOffset = pset.get<size_t >("CoherentOffset", 0);
fStructuringElement = pset.get<size_t >("StructuringElement", 20);
fMorphWindow = pset.get<size_t >("FilterWindow", 10);
Expand Down Expand Up @@ -242,7 +241,8 @@ void TPCNoiseFilter1DMC::configure(fhicl::ParameterSet const &pset)

void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&,
const daq::INoiseFilter::ChannelPlaneVec& channelPlaneVec,
const icarus_signal_processing::ArrayFloat& dataArray)
const icarus_signal_processing::ArrayFloat& dataArray,
const size_t& coherentNoiseGrouping)
{
cet::cpu_timer theClockTotal;

Expand All @@ -268,7 +268,7 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&,
if (fNumTruncBins.size() < numChannels) fNumTruncBins.resize(numChannels);
if (fRangeBins.size() < numChannels) fRangeBins.resize(numChannels);

if (fThresholdVec.size() < numChannels) fThresholdVec.resize(numChannels / fCoherentNoiseGrouping);
if (fThresholdVec.size() < numChannels) fThresholdVec.resize(numChannels / coherentNoiseGrouping);

if (fFilterFunctionVec.size() < numChannels) fFilterFunctionVec.resize(numChannels);

Expand All @@ -287,10 +287,11 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&,
std::vector<geo::WireID> widVec = fGeometry->ChannelToWire(fChannelIDVec[idx]);

// Handle the filter function to use for this channel
unsigned int plane = channelPlaneVec[idx].second;
// Note the modulus... this to enable a workaround for the wirecell 2D drift which misses the channels with no signal
unsigned int plane = channelPlaneVec[idx].second % 3;

// Set the threshold which toggles between planes
fThresholdVec[idx / fCoherentNoiseGrouping] = fThreshold[plane];
fThresholdVec[idx / coherentNoiseGrouping] = fThreshold[plane];

switch(fFilterModeVec[plane][0])
{
Expand Down Expand Up @@ -341,7 +342,7 @@ void TPCNoiseFilter1DMC::process_fragment(detinfo::DetectorClocksData const&,
fFilterFunctionVec.begin(),
fThresholdVec,
numChannels,
fCoherentNoiseGrouping,
coherentNoiseGrouping,
fCoherentNoiseOffset,
fMorphWindow);

Expand Down
13 changes: 7 additions & 6 deletions icaruscode/Decode/DecoderTools/TPCNoiseFilterCanny_tool.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,8 @@ class TPCNoiseFilterCannyMC : virtual public INoiseFilter
*/
virtual void process_fragment(detinfo::DetectorClocksData const&,
const daq::INoiseFilter::ChannelPlaneVec&,
const icarus_signal_processing::ArrayFloat&) override;
const icarus_signal_processing::ArrayFloat&,
const size_t&) override;

/**
* @brief Recover the channels for the processed fragment
Expand Down Expand Up @@ -274,7 +275,6 @@ void TPCNoiseFilterCannyMC::configure(fhicl::ParameterSet const &pset)

fMorphologicalFilter = std::make_unique<icarus_signal_processing::Dilation2D>(fMorph2DStructuringElementX,fMorph2DStructuringElementY);

fCoherentNoiseGrouping = pset.get<unsigned int >("CoherentNoiseGrouping", 32);
fCoherentNoiseOffset = pset.get<unsigned int >("CoherentNoiseOffset", 24);
fMorphologicalWindow = pset.get<unsigned int >("MorphologicalWindow", 10);
fCoherentThresholdFactor = pset.get<float >("CoherentThresholdFactor", 2.5);
Expand Down Expand Up @@ -330,8 +330,9 @@ void TPCNoiseFilterCannyMC::configure(fhicl::ParameterSet const &pset)
}

void TPCNoiseFilterCannyMC::process_fragment(detinfo::DetectorClocksData const&,
const daq::INoiseFilter::ChannelPlaneVec& channelPlaneVec,
const icarus_signal_processing::ArrayFloat& dataArray)
const daq::INoiseFilter::ChannelPlaneVec& channelPlaneVec,
const icarus_signal_processing::ArrayFloat& dataArray,
const size_t& coherentNoiseGrouping)
{
cet::cpu_timer theClockTotal;

Expand All @@ -357,7 +358,7 @@ void TPCNoiseFilterCannyMC::process_fragment(detinfo::DetectorClocksData const&,
if (fNumTruncBins.size() < numChannels) fNumTruncBins.resize(numChannels);
if (fRangeBins.size() < numChannels) fRangeBins.resize(numChannels);

if (fThresholdVec.size() < numChannels) fThresholdVec.resize(numChannels / fCoherentNoiseGrouping);
if (fThresholdVec.size() < numChannels) fThresholdVec.resize(numChannels / coherentNoiseGrouping);

if (fFilterFunctionVec.size() < numChannels) fFilterFunctionVec.resize(numChannels);

Expand All @@ -381,7 +382,7 @@ void TPCNoiseFilterCannyMC::process_fragment(detinfo::DetectorClocksData const&,
unsigned int plane = widVec[0].Plane;

// Set the threshold which toggles between planes
fThresholdVec[idx / fCoherentNoiseGrouping] = fThreshold[plane];
fThresholdVec[idx / coherentNoiseGrouping] = fThreshold[plane];

switch(fFilterModeVec[plane][0])
{
Expand Down
1 change: 0 additions & 1 deletion icaruscode/Decode/DecoderTools/decoderTools_icarus.fcl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,6 @@ TPCNoiseFilter1DTool: {
tool_type: TPCNoiseFilter1D
fragment_id_offset: 0
NSigmaForTrucation: 3.5
CoherentGrouping: 64
StructuringElement: 16
FilterWindow: 10
Threshold: [12.0, 12.0, 12.0] #[8., 7.0 , 8.0] #[2.75, 2.75, 2.75] # changing threshold scheme from dynamic to static
Expand Down
Loading

0 comments on commit 4403b63

Please sign in to comment.